simple geometric model

relative invagination height

\begin{equation} \label{rel-h} h = \frac{H_{i}}{\left(\frac{S_{i}}{\pi}\right)^{1/2}} \end{equation} \begin{equation} H_{i} = h \left(\frac{S_{i}}{\pi}\right)^{1/2} \end{equation}

droplet-membrane contact angle

\begin{equation} \begin{array}{l} \gamma_{S G}=\gamma_{S L}+\gamma_{L G} \cos \theta \\ \sigma+\gamma_{m}=\sigma+\gamma_{d} \cos \theta_{d} \\ \cos \theta_{d}=\frac{\gamma_{m}}{\gamma_{d}} \\ \sin \theta_{d}=\sqrt{1 - \cos \theta_{d}^{2}} \\ \sin \theta_{d}=\sqrt{1 - (\frac{\gamma_{m}}{\gamma_{d}})^{2}} \end{array} \end{equation} \begin{equation} \label{droplet-contact-angle} \theta_{d}=\arccos \left(\frac{\gamma_{m}}{\gamma_{d}}\right) \end{equation}

outside membrane area

circular plane radii

droplet-cytosol contact area

invagination surface area

invagination volume

invagination plane radius

invagination radius

invagination contact angle

TODO droplet volume

droplet + invagination volume

droplet + invagination spherical cap radius

droplet + invagination height

flat membrane-droplet contact area

full derivations

• \(F_{1}(\Delta S_{o})\)
  
  \begin{equation} F_{6}(S_{o}) = \sigma S_{o} \end{equation}

  substituting outside-membrane-area:

  \begin{equation} F_{6}(S_{o}) = \sigma \pi\left(A_{t}^{2}-A_{d}^{2}\right) \end{equation}

  substituting droplet-plane-radius:

  \begin{equation} F_{6}(S_{o}) = \sigma \pi\left(A_{t}^{2}-(R_{d} \sin \theta_{d})^{2}\right) \end{equation}

  substituting di-cap-radius:

  \begin{equation} F_{6}(S_{o}) = \sigma \pi\left(A_{t}^{2}-\left(\left(\frac{3 V_{b}}{\pi\left(2+\cos \theta_{d}\right)\left(1-\cos \theta_{d}\right)^{2}}\right)^{1 / 3} \sin \theta_{d}\right)^{2}\right) \end{equation}

  substituting droplet-contact-angle:

  \begin{equation} F_{6}(S_{o}) = \sigma \pi\left(A_{t}^{2}-\left(\left(\frac{3 V_{b}}{\pi\left(2+\frac{\gamma_{m}}{\gamma_{d}}\right)\left(1-\frac{\gamma_{m}}{\gamma_{d}}\right)^{2}}\right)^{1 / 3} \sqrt{1 - (\frac{\gamma_{m}}{\gamma_{d}})^{2}}\right)^{2}\right) \end{equation}

  substituting di-volume:

  \begin{equation} F_{6}(S_{o}) = \sigma \pi\left(A_{t}^{2}-\left(\left(\frac{3 (V_{d}+V_{i})}{\pi\left(2+\frac{\gamma_{m}}{\gamma_{d}}\right)\left(1-\frac{\gamma_{m}}{\gamma_{d}}\right)^{2}}\right)^{1 / 3} \sqrt{1 - (\frac{\gamma_{m}}{\gamma_{d}})^{2}}\right)^{2}\right) \end{equation}

  substituting invagination-volume:

  \begin{equation} F_{6}(S_{o}) = \sigma \pi\left(A_{t}^{2}-\left(\left(\frac{3 \left(V_{d}+\frac{\pi}{3} H_{i}^{2} \left(3 R_{i} - H_{i} \right)\right)}{\pi\left(2+\frac{\gamma_{m}}{\gamma_{d}}\right)\left(1-\frac{\gamma_{m}}{\gamma_{d}}\right)^{2}}\right)^{1 / 3} \sqrt{1 - (\frac{\gamma_{m}}{\gamma_{d}})^{2}}\right)^{2}\right) \end{equation}

  substituting invagination-radius:

  \begin{equation} F_{6}(S_{o}) = \sigma \pi\left(A_{t}^{2}-\left(\left(\frac{3 \left(V_{d}+\frac{\pi}{3} H_{i}^{2} \left(3 \frac{S_{i}}{2 \pi H_{i}} - H_{i} \right)\right)}{\pi\left(2+\frac{\gamma_{m}}{\gamma_{d}}\right)\left(1-\frac{\gamma_{m}}{\gamma_{d}}\right)^{2}}\right)^{1 / 3} \sqrt{1 - (\frac{\gamma_{m}}{\gamma_{d}})^{2}}\right)^{2}\right) \end{equation}

  “simplified” in mathematica:

  \begin{equation} F_{6}(S_{o}) = \sigma \pi \left(A_t^2+\frac{\left(\gamma _m^2-\gamma _d^2\right) \left(\frac{\gamma _d^2 \left(3 V_d+\pi H_i^2\left(\frac{3 S_i}{2 \pi H_i}-H_i\right)\right)}{\pi \left(\frac{\gamma _m}{\gamma _d}+2\right) \left(\gamma _d-\gamma _m\right)^{2}}\right)^{2/3}}{\gamma _d^2}\right) \end{equation}

  rearranged to isolate constants vs. functions of \(H_{i}\):

  

  \begin{equation} F_{6}(S_{o}) = \sigma S_{o} \end{equation} \begin{equation} S_{o} = \Delta S_{o} + \pi A_{t}^{2} \end{equation} \begin{equation} F_{6}(S_{o}) = F_{1}(\Delta S_{o}) + \sigma \pi A_{t}^{2} \end{equation} \begin{equation} F_{1}(\Delta S_{o}) = \sigma \Delta S_{o} \end{equation} \begin{equation} \Delta S_{o} = S_{o} - \pi A_{t}^{2} \end{equation} \begin{equation} \Delta S_{o}=\left(\frac{\pi\left(\gamma_{m}^{2}-\gamma_{d}^{2}\right)^{3}}{\left(2 \gamma_{d}+\gamma_{m}\right)^{2}\left(\gamma_{d}-\gamma_{m}\right)^{4}}\right)^{1 / 3}\left(-\pi H_{i}^{3}+\frac{3}{2} S_{i} H_{i}+3 V_{d}\right)^{2 / 3} \end{equation}

  replacing invagination height from rel-h:

  \begin{equation} H_{i} = h \left(\frac{S_{i}}{\pi}\right)^{1/2} \end{equation} \begin{equation} \Delta S_{o}=\left(\frac{\pi\left(\gamma_{m}^{2}-\gamma_{d}^{2}\right)^{3}}{\left(2 \gamma_{d}+\gamma_{m}\right)^{2}\left(\gamma_{d}-\gamma_{m}\right)^{4}}\right)^{1 / 3}\left(-\pi \left(h \left(\frac{S_{i}}{\pi}\right)^{1/2}\right)^{3}+\frac{3}{2} S_{i} h \left(\frac{S_{i}}{\pi}\right)^{1/2}+3 V_{d}\right)^{2 / 3} \end{equation}

  simplified:

  \begin{equation} \Delta S_{o}=\left(\frac{\pi\left(\gamma_{m}^{2}-\gamma_{d}^{2}\right)^{3}}{\left(2 \gamma_{d}+\gamma_{m}\right)^{2}\left(\gamma_{d}-\gamma_{m}\right)^{4}}\right)^{1 / 3}\left(\frac{S_{i}^{3 / 2}}{\sqrt{\pi}}\left(\frac{3}{2} h-h^{3}\right)+3 V_{d}\right)^{2/3} \end{equation} \begin{equation} \Delta S_{o}=\alpha \left(\frac{S_{i}^{3 / 2}}{\sqrt{\pi}}\left(\frac{3}{2} h-h^{3}\right)+3 V_{d}\right)^{2/3} \end{equation} \begin{equation} \alpha = \left(\frac{\pi\left(\gamma_{m}^{2}-\gamma_{d}^{2}\right)^{3}}{\left(2 \gamma_{d}+\gamma_{m}\right)^{2}\left(\gamma_{d}-\gamma_{m}\right)^{4}}\right)^{1 / 3} \end{equation}
• \(F_{2}(\Delta S_{m})\)
  

  gathering from free-energy:

  \begin{equation} \label{} F_{7}(S_{m}) = \left(\sigma - \gamma_{m} \right) S_{m} \end{equation}

  from md-area:

  \begin{equation} \label{md-area} S_{m}=\pi\left(A_{d}^{2}-A_{i}^{2}\right) \end{equation}

  substitute droplet-plane-radius:

  \begin{equation} A_{d}=R_{d} \sin \theta_{d} \end{equation}

  and invagination-plane-radius:

  \begin{equation} A_{i}=\sqrt{\frac{S_{i}}{\pi}-H_{i}^{2}} \end{equation}

  to get:

  \begin{equation} S_{m}=\pi\left(\left(R_{d} \sin \theta_{d}\right)^{2}-\frac{S_{i}}{\pi}+H_{i}^{2}\right) \end{equation}

  rearranged:

  \begin{equation} S_{m}=\pi\left(\left(R_{d} \sin \theta_{d}\right)^{2}+H_{i}^{2}\right) - S_{i} \end{equation}

  substitute di-cap-radius:

  \begin{equation} R_{d}=\left(\frac{3 V_{b}}{\pi\left(2+\cos \theta_{d}\right)\left(1-\cos \theta_{d}\right)^{2}}\right)^{1 / 3} \end{equation}

  to get:

  \begin{equation} S_{m}=\pi\left(\left(\left(\frac{3 V_{b}}{\pi\left(2+\cos \theta_{d}\right)\left(1-\cos \theta_{d}\right)^{2}}\right)^{1 / 3} \sin \theta_{d}\right)^{2}+H_{i}^{2}\right) - S_{i} \end{equation}

  substitute di-volume:

  \begin{equation} V_{b}=V_{d}+V_{i} \end{equation}

  to get:

  \begin{equation} S_{m}=\pi\left(\left(\left(\frac{3 \left(V_{d}+V_{i} \right)}{\pi\left(2+\cos \theta_{d}\right)\left(1-\cos \theta_{d}\right)^{2}}\right)^{1 / 3} \sin \theta_{d}\right)^{2}+H_{i}^{2}\right) - S_{i} \end{equation}

  substitute invagination-volume:

  \begin{equation} V_{i}=\frac{\pi}{3} H_{i}^{2} \left(3 R_{i} - H_{i} \right) \end{equation}

  to get:

  \begin{equation} S_{m}=\pi\left(\left(\left(\frac{3 \left(V_{d}+\frac{\pi}{3} H_{i}^{2} \left(3 R_{i} - H_{i} \right) \right)}{\pi\left(2+\cos \theta_{d}\right)\left(1-\cos \theta_{d}\right)^{2}}\right)^{1 / 3} \sin \theta_{d}\right)^{2}+H_{i}^{2}\right) - S_{i} \end{equation}

  substitute invagination-radius:

  \begin{equation} R_{i} =\frac{S_{i}}{2 \pi H_{i}} \end{equation}

  to get:

  \begin{equation} S_{m}=\pi\left(\left(\left(\frac{3 \left(V_{d}+\frac{\pi}{3} H_{i}^{2} \left(3 \frac{S_{i}}{2 \pi H_{i}} - H_{i} \right) \right)}{\pi\left(2+\cos \theta_{d}\right)\left(1-\cos \theta_{d}\right)^{2}}\right)^{1 / 3} \sin \theta_{d}\right)^{2}+H_{i}^{2}\right) - S_{i} \end{equation}

  substitute droplet-contact-angle:

  \begin{equation} \begin{array}{l} \cos \theta_{d}=\frac{\gamma_{m}}{\gamma_{d}} \\ \sin \theta_{d}=\sqrt{1 - (\frac{\gamma_{m}}{\gamma_{d}})^{2}} \end{array} \end{equation}

  to get:

  \begin{equation} S_{m}=\pi\left(\left(\left(\frac{3 \left(V_{d}+\frac{\pi}{3} H_{i}^{2} \left(3 \frac{S_{i}}{2 \pi H_{i}} - H_{i} \right) \right)}{\pi\left(2+\frac{\gamma_{m}}{\gamma_{d}}\right)\left(1-\frac{\gamma_{m}}{\gamma_{d}}\right)^{2}}\right)^{1 / 3} \sqrt{1 - (\frac{\gamma_{m}}{\gamma_{d}})^{2}}\right)^{2}+H_{i}^{2}\right) - S_{i} \end{equation}

  simplified in mathematica:

  \begin{equation} S_{m}=\pi \left(\frac{\left(\gamma _d^2-\gamma _m^2\right) \left(\frac{\gamma _d^2 \left(3 V_d+\pi H_i^2\left(\frac{3 S_i}{2 \pi H_i}-H_i\right)\right)}{\pi \left(2-\frac{\gamma _m}{\gamma _d}\right) \left(\gamma _d+\gamma _m\right)^{2}}\right)^{2/3}}{\gamma _d^2}+H_i^2\right)-S_i \end{equation}

  manually rearranged:

  

  \begin{equation} \label{} F_{7}(S_{m}) = \left(\sigma - \gamma_{m} \right) S_{m} \end{equation} \begin{equation} S_{m}= \Delta S_{m} - S_{i} \end{equation} \begin{equation} \label{} F_{7}(S_{m}) = F_{2}(\Delta S_{m}) - \left(\sigma - \gamma_{m} \right) S_{i} \end{equation} \begin{equation} \label{} F_{2}(\Delta S_{m}) = \left(\sigma - \gamma_{m} \right) \Delta S_{m} \end{equation} \begin{equation} \Delta S_{m}=\left(\frac{\pi\left(\gamma_{d}^{2}-\gamma_{m}^{2}\right)^{3}}{\left(2 \gamma_{d}+\gamma_{m}\right)^{2}\left(\gamma_{d}-\gamma_{m}\right)^{4}}\right)^{1 / 3}\left(-\pi H_{i}^{3}+\frac{3}{2} S_{i} H_{i}+3 V_{d}\right)^{2 / 3}+\pi H_{i}^{2} \end{equation}

  replacing invagination height from rel-h:

  \begin{equation} H_{i} = h \left(\frac{S_{i}}{\pi}\right)^{1/2} \end{equation} \begin{equation} \Delta S_{m}=\left(\frac{\pi\left(\gamma_{d}^{2}-\gamma_{m}^{2}\right)^{3}}{\left(2 \gamma_{d}+\gamma_{m}\right)^{2}\left(\gamma_{d}-\gamma_{m}\right)^{4}}\right)^{1 / 3}\left(-\pi \left(h \left(\frac{S_{i}}{\pi}\right)^{1/2}\right)^{3}+\frac{3}{2} S_{i} \left(h \left(\frac{S_{i}}{\pi}\right)^{1/2}\right)+3 V_{d}\right)^{2 / 3}+ \pi \left(h \left(\frac{S_{i}}{\pi}\right)^{1/2}\right)^{2} \end{equation}

  simplified:

  \begin{equation} \Delta S_{m}=\left(\frac{\pi\left(\gamma_{d}^{2}-\gamma_{m}^{2}\right)^{3}}{\left(2 \gamma_{d}+\gamma_{m}\right)^{2}\left(\gamma_{d}-\gamma_{m}\right)^{4}}\right)^{1 / 3} \left(\frac{S_{i}^{3 / 2}}{\sqrt{\pi}}\left(\frac{3}{2} h-h^{3}\right)+3 V_{d}\right)^{2/3} + S_{i} h^{2} \end{equation} \begin{equation} \Delta S_{m}= \beta \left(\frac{S_{i}^{3 / 2}}{\sqrt{\pi}}\left(\frac{3}{2} h-h^{3}\right)+3 V_{d}\right)^{2/3} + S_{i} h^{2} \end{equation} \begin{equation} \beta = \left(\frac{\pi\left(\gamma_{d}^{2}-\gamma_{m}^{2}\right)^{3}}{\left(2 \gamma_{d}+\gamma_{m}\right)^{2}\left(\gamma_{d}-\gamma_{m}\right)^{4}}\right)^{1 / 3} \end{equation}
• \(F_{3}(S_{d})\)
  

  gathering from free-energy:

  \begin{equation} \label{} F_{3}(S_{d}) = \gamma_{d} S_{d} \end{equation}

  from dc-contact-area:

  \begin{equation} S_{d}=2 \pi R_{d}^{2}\left(1-\cos \theta_{d}\right) \end{equation}

  substitute di-cap-radius

  \begin{equation} R_{d}=\left(\frac{3 V_{b}}{\pi\left(2+\cos \theta_{d}\right)\left(1-\cos \theta_{d}\right)^{2}}\right)^{1 / 3} \end{equation}

  to get:

  \begin{equation} S_{d}=2 \pi \left(\frac{3 V_{b}}{\pi\left(2+\cos \theta_{d}\right)\left(1-\cos \theta_{d}\right)^{2}}\right)^{2 / 3}\left(1-\cos \theta_{d}\right) \end{equation}

  substitute di-volume

  \begin{equation} V_{b}=V_{d}+V_{i} \end{equation}

  to get:

  \begin{equation} S_{d}=2 \pi \left(\frac{3 \left(V_{d}+V_{i}\right)}{\pi\left(2+\cos \theta_{d}\right)\left(1-\cos \theta_{d}\right)^{2}}\right)^{2 / 3}\left(1-\cos \theta_{d}\right) \end{equation}

  substitute invagination-volume

  \begin{equation} V_{i}=\frac{\pi}{3} H_{i}^{2} \left(3 R_{i} - H_{i} \right) \end{equation}

  to get:

  \begin{equation} S_{d}=2 \pi \left(\frac{3 \left(V_{d}+\frac{\pi}{3} H_{i}^{2} \left(3 R_{i} - H_{i} \right)\right)}{\pi\left(2+\cos \theta_{d}\right)\left(1-\cos \theta_{d}\right)^{2}}\right)^{2 / 3}\left(1-\cos \theta_{d}\right) \end{equation}

  substitute droplet-contact-angle

  \begin{equation} \cos \theta_{d}=\frac{\gamma_{m}}{\gamma_{d}} \end{equation}

  to get:

  \begin{equation} S_{d}=2 \pi \left(\frac{3 \left(V_{d}+\frac{\pi}{3} H_{i}^{2} \left(3 R_{i} - H_{i} \right)\right)}{\pi\left(2+\frac{\gamma_{m}}{\gamma_{d}}\right)\left(1-\frac{\gamma_{m}}{\gamma_{d}}\right)^{2}}\right)^{2 / 3}\left(1-\frac{\gamma_{m}}{\gamma_{d}}\right) \end{equation}

  substitute invagination-radius

  \begin{equation} R_{i} &=\frac{S_{i}}{2 \pi H_{i}} \end{equation}

  to get:

  \begin{equation} S_{d}=2 \pi \left(\frac{3 \left(V_{d}+\frac{\pi}{3} H_{i}^{2} \left(3 \frac{S_{i}}{2 \pi H_{i}} - H_{i} \right)\right)}{\pi\left(2+\frac{\gamma_{m}}{\gamma_{d}}\right)\left(1-\frac{\gamma_{m}}{\gamma_{d}}\right)^{2}}\right)^{2 / 3}\left(1-\frac{\gamma_{m}}{\gamma_{d}}\right) \end{equation}

  simplified in mathematica:

  \begin{equation} S_{d}= \frac{2 \pi \left(\gamma _d-\gamma _m\right) \left(\frac{\gamma _d^2 \left(3 V_d+\pi H_i^2\left(\frac{3 S_i}{2 \pi H_i}-H_i\right)\right)}{\pi \left(\frac{\gamma _m}{\gamma _d}+2\right) \left(\gamma _d-\gamma _m\right)^{2}}\right)^{2/3}}{\gamma _d} \end{equation}

  manually rearranged:

  

  \begin{equation} S_{d}=\left(\frac{8 \pi}{\left(2+\frac{\gamma_{m}}{\gamma_{d}}\right)^{2}\left(1-\frac{\gamma_{m}}{\gamma_{d}}\right)}\right)^{1 / 3}\left(\frac{S_{i}^{3 / 2}}{\sqrt{\pi}}\left(\frac{3}{2} h-h^{3}\right)+3 V_{d}\right)^{2/3} \end{equation} \begin{equation} S_{d}=\epsilon \left(\frac{S_{i}^{3 / 2}}{\sqrt{\pi}}\left(\frac{3}{2} h-h^{3}\right)+3 V_{d}\right)^{2/3} \end{equation} \begin{equation} \epsilon = \left(\frac{8 \pi}{\left(2+\frac{\gamma_{m}}{\gamma_{d}}\right)^{2}\left(1-\frac{\gamma_{d}}{\gamma_{m}}\right)}\right)^{1 / 3} \end{equation}
• \(F_{4}(R_{i})\)
  

  gathering from f-forces:

  \begin{equation} \label{} \begin{array}{l} F_{4}(R_{i}) = \frac{K}{2} S_{i} \left(\frac{2}{R_{i}} - 2 C_{i}\right)^{2} \end{array} \end{equation}

  from invagination-radius

  \begin{equation} R_{i} &=\frac{S_{i}}{2 \pi H_{i}} \end{equation}

  replacing invagination height from rel-h:

  \begin{equation} H_{i} = h \left(\frac{S_{i}}{\pi}\right)^{1/2} \end{equation} \begin{equation} R_{i} = \frac{S_{i}}{2 \pi h \left(\frac{S_{i}}{\pi}\right)^{1/2}} \end{equation}

  simplified:

  \begin{equation} R_{i} = \frac{\sqrt{S_i}}{2 \sqrt{\pi } h} \end{equation}
  • \(F_{4}(h)\)
    

    substitute invagination-radius

    \begin{equation} R_{i} &=\frac{S_{i}}{2 \pi H_{i}} \end{equation}

    to get:

    \begin{equation} \label{} \begin{array}{l} F_{4}(H_{i}) = \frac{K}{2} S_{i} \left(\frac{2}{\frac{S_{i}}{2 \pi H_{i}}} - 2 C_{i}\right)^{2} \end{array} \end{equation}

    simplified:

    \begin{equation} \label{} F_{4}(H_{i}) = \frac{2 K}{S_i} \left(C_i S_i-2 \pi H_i\right)^{2} \end{equation}

    replacing invagination height from rel-h:

    \begin{equation} H_{i} = h \left(\frac{S_{i}}{\pi}\right)^{1/2} \end{equation} \begin{equation} \label{} F_{4}(h) = \frac{2 K}{S_i} \left(C_i S_i-2 \pi \left(h \left(\frac{S_{i}}{\pi}\right)^{1/2}\right)\right)^{2} \end{equation}

    simplified:

    \begin{equation} \label{} F_{4}(h) = 2 K \left(C_i \sqrt{S_i}-2 \sqrt{\pi } h\right)^2 \end{equation}
• \(F_{5}(V_{i})\)
  

  gathering from f-forces:

  \begin{equation} \label{} F_{5}(V_{i}) = P V_{i} \end{equation}

  from invagination-volume-a-h:

  \begin{equation} \label{} \begin{array}{l} V_{i}=\frac{1}{6} \pi H_{i}\left(3 A_{i}^{2}+H_{i}^{2}\right) \end{array} \end{equation}

  substitute invagination-plane-radius

  \begin{equation} \label{} \begin{array}{l} A_{i}=\sqrt{\frac{S_{i}}{\pi}-H_{i}^{2}} \end{array} \end{equation}

  to get:

  \begin{equation} \label{} \begin{array}{l} V_{i}=\frac{1}{6} \pi H_{i}\left(3 \left(\frac{S_{i}}{\pi}-H_{i}^{2}\right)+H_{i}^{2}\right) \end{array} \end{equation}

  simplified:

  \begin{equation} \label{} \begin{array}{l} V_{i}=\frac{1}{6} \pi H_{i}\left(\frac{3 S_i}{\pi }-2 H_{i}^2\right) \end{array} \end{equation}

  replacing invagination height from rel-h:

  \begin{equation} H_{i} = h \left(\frac{S_{i}}{\pi}\right)^{1/2} \end{equation} \begin{equation} \label{} V_{i}=\frac{1}{6} \pi h \left(\frac{S_{i}}{\pi}\right)^{1/2}\left(\frac{3 S_i}{\pi }-2 \left(h \left(\frac{S_{i}}{\pi}\right)^{1/2}\right)^2\right) \end{equation}

  simplified:

  \begin{equation} \label{} V_{i}= \frac{h \left(3-2 h^2\right) S_i^{3/2}}{6 \sqrt{\pi }} \end{equation} \begin{equation} \label{} V_{i}=\eta \left(3 h-2 h^{3}\right) \end{equation} \begin{equation} \label{} \eta=\frac{S_i^{3/2}}{6 \sqrt{\pi }} \end{equation}
• \(D\): everything else that doesn’t change with \(H_{i}\)
  
  \begin{equation} \label{D} D = \sigma S_{i} - \gamma_{i} S_{i} \end{equation}
• \(C\): additional constants grouped out from \(\Delta\) terms
  
  \begin{equation} \label{} C = D + \sigma \pi A_{t}^{2} - \left(\sigma - \gamma_{m} \right) S_{i} \end{equation} \begin{equation} \label{} D = \sigma S_{i} - \gamma_{i} S_{i} \end{equation}

  simplified:

  \begin{equation} \label{} C = \sigma \pi A_t^2+S_i \left(\gamma _m-\gamma _i\right) \end{equation}

DONE substituted after isolating \(h\)

• could anything be factored out more first?

\begin{equation} \label{} F(h) = \sigma \alpha \xi + \left(\sigma - \gamma_{m} \right) \left( \beta \xi + S_{i} h^{2}\right) + \gamma_{d} \epsilon \xi + 2 K \left(C_i \sqrt{S_i}-2 \sqrt{\pi } h\right)^{2} + P \eta \left(3 h-2 h^{3}\right) + C \end{equation} \begin{equation} \label{} \xi = \left(\frac{S_{i}^{3 / 2}}{\sqrt{\pi}}\left(\frac{3}{2} h-h^{3}\right)+3 V_{d}\right)^{2/3} \end{equation}

\begin{equation} F(h)=\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \epsilon\right)\left(\frac{S_{i}^{3 / 2}}{\sqrt{\pi}}\left(\frac{3}{2} h-h^{3}\right)+3 V_{d}\right)^{2/3}+\left(\sigma-\gamma_{m}\right) S_{i} h^{2} +2 K\left(C_{i} \sqrt{S}_{i}-2 \sqrt{\pi} h\right)^{2}+P \eta \left(3 h-2 h^{3}\right)+C \end{equation} \begin{equation} F(h)= \left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \epsilon\right) \left(\frac{S_{i}^{3 / 2}}{\sqrt{\pi}}\left(\frac{3}{2} h-h^{3}\right)+3 V_{d}\right)^{2/3} - 2 P \eta h^{3}+\left(\left(\sigma-\gamma_{m}\right) S_{i}+8 K \pi\right) h^{2} +\left(3 P \eta-8 C_{i} \sqrt{\pi S_{i}}\right) h+B \end{equation} \begin{equation} B=C+2 K C_{i}^{2} S_{i} \end{equation}

testing sign of alpha terms

import numpy as np
import math
def cube_root(i):
    return math.pow(abs(i),float(1)/3) * (1,-1)[i<0]

def alpha_so(gamma_m, gamma_d):
    alpha = cube_root((np.pi*(gamma_m**2 - gamma_d**2)**3)/((2*gamma_d+gamma_m)**2 * (gamma_d-gamma_m)**4))
    return alpha

for gamma_m in np.linspace(0, 200, 9):
    for gamma_d in np.linspace(0, 200, 9):
        print('gamma_m',gamma_m,'gamma_d',gamma_d,alpha_so(gamma_m, gamma_d))

thermodynamic spontaneity (global stability) [5/5]

• DONE setting up equality
  

  

  \begin{equation} \left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \in\right)\left(3 V_{d}\right)^{2 / 3}+B=\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \in\right)\left(3 V_d+\frac{S_i^{3/2}}{2 \sqrt{\pi }}\right)^{2/3} + \left(\sigma-\gamma_{m}\right) S_{i}+8 K \pi+ P \eta-8 K C_{i} \sqrt{\pi S_{i}}+B \end{equation}
• DONE solving for \(K\)
  

  

  \begin{equation} \label{therm_spont_k} K = \frac{\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \in\right)\left(\left(3 V_{d}\right)^{2 / 3}-\left(\frac{S_{i}^{3 / 2}}{2 \sqrt{\pi}}+3 V_{d}\right)^{2 / 3}\right)-\left(\sigma-\gamma_{m}\right) S_{i}- P \eta}{8 \pi-8 C_{i} \sqrt{\pi S_{i}}} \end{equation}
• DONE solving for \(P\)
  

  

  \begin{equation} \label{therm_spont_p} P=\frac{\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \in\right)\left(\left(3 V_{d}\right)^{2 / 3}-\left(\frac{S_{i}^{3 / 2}}{2 \sqrt{\pi}}+3 V_{d}\right)^{2 / 3}\right)-\left(\sigma-\gamma_{m}\right) S_{i}-8 K \left(\pi - C_{i} \sqrt{\pi S_{i}}\right)}{\eta} \end{equation}
• DONE solving for \(\sigma\)
  

  

  \begin{equation} \label{therm_spont_sigma} \sigma=\frac{\frac{-\gamma_{m} S_{i}+8 K \left(\pi - C_{i} \sqrt{\pi S_{i}}\right) + P \eta}{\left(3 V_{d}\right)^{2 / 3}-\left(\frac{S_{i}^{3 / 2}}{2 \sqrt{\pi}}+3 V_{d}\right)^{2 / 3}}-\gamma_{d} \epsilon+\gamma_{m} \beta}{\alpha+\beta-\frac{S_{i}}{\left(3 V_{d}\right)^{2 / 3}-\left(\frac{S_{i}^{3 / 2}}{2 \sqrt{\pi}}+3 V_{d}\right)^{2 / 3}}} \end{equation}
• DONE solving for \(C_{i}\)
  

  

  \begin{equation} \label{therm_spont_ci} C_{i} = \frac{\pi-\frac{\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \in\right)\left(\left(3 V_{d}\right)^{2 / 3}-\left(\frac{S_{i}^{3 / 2}}{2 \sqrt{\pi}}+3 V_{d}\right)^{2 / 3}\right)-\left(\sigma-\gamma_{m}\right) S_{i}-P_{\eta}}{8 K}}{\sqrt{\pi S_{i}}} \end{equation}
• any others?
  

kinetic spontaneity (local stability) [11/11]

• DONE solving derivative
  

  

  \begin{equation} \label{f-prime} F^{\prime}(h)=\frac{\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \epsilon\right)\left(1-2 h^{2}\right) S_{i}^{3 / 2}}{\sqrt{\pi}\left(\frac{\left(\frac{3 h}{2}-h^{3}\right) S_{i}^{3 / 2}}{\sqrt{\pi}}+3 V_{d}\right)^{1 / 3}}-6 P \eta h^{2}+2\left(\left(\sigma-\gamma_{m}\right) S_{i}+8 K \pi\right) h+3 P \eta-8 K C_{i} \sqrt{\pi S_{i}} \end{equation}
• DONE solving for h=0
  
  \begin{equation} \label{f-prime-alpha-beta} F^{\prime}(0)=\frac{\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \epsilon\right) S_{i}^{3 / 2}}{\sqrt{\pi} \left(3 V_{d}\right)^{1 / 3}}+3 P \eta-8 K C_{i} \sqrt{\pi S_{i}} \end{equation} \begin{equation} \label{f-prime} F^{\prime}(0)=\frac{\left(\gamma_{m} \alpha + \gamma_{d} \epsilon \right) S_{i}^{3 / 2}}{\sqrt{\pi} \left(3 V_{d}\right)^{1 / 3}}+3 P \eta-8 K C_{i} \sqrt{\pi S_{i}} \end{equation}
• DONE sign of droplet prefactor at h=0
  

  

  the sign of the term concerning the droplet must always be positive.

• DONE solving for \(K\)
  

  

  \begin{equation} \label{kinetic_spont_k} K=\frac{\frac{\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \epsilon\right) S_{i}^{3 / 2}}{\sqrt{\pi}\left(3 V_{d}\right)^{1 / 3}}+3 P \eta}{8 C_{i} \sqrt{\pi S_{i}}} \end{equation}
• DONE solving for \(P\)
  

  

  \begin{equation} \label{kinetic_spont_p} P=\frac{\frac{-\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \epsilon\right) S_{i}^{3 / 2}}{\sqrt{\pi}\left(3 V_{d}\right)^{1 / 3}}+8 K C_{i} \sqrt{\pi S_{i}}}{3 \eta} \end{equation}
• DONE solving for \(\sigma\)
  

  

  \begin{equation} \label{kinetic_spont_sigma} \sigma=\frac{\frac{\sqrt{\pi}\left(3 V_{d}\right)^{1 / 3}}{S_{i}^{3 / 2}}\left(8 K C_{i} \sqrt{\pi S_{i}}-3 P \eta \right)-\gamma_{d} \epsilon+\gamma_{m} \beta}{(\alpha+\beta)} \end{equation}
• DONE solving for \(C_{i}\)
  

  

  \begin{equation} \label{kinetic_spont_ci} C_{i} = \frac{\frac{\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \epsilon\right) S_{i}^{3 / 2}} {\sqrt{\pi}\left(3 V_{d}\right)^{1 / 3}}+3 P \eta}{8 K \sqrt{\pi S_{i}}} \end{equation}
• CANCELLED solving for \(\gamma_{m}\)
  

  

  \begin{equation} \label{kinetic_spont_gammam} \gamma_{m}=\frac{\frac{-\sqrt{\pi}\left(3 V_{d}\right)^{1 / 3}}{S_{i}^{3 / 2}}\left(8 K C_{i} \sqrt{\pi S_{i}}-3 P \eta \right)+\gamma_{d} \epsilon+\sigma \alpha+\sigma \beta}{\beta} \end{equation}
• CANCELLED solving for \(\gamma_{d}\)
  

  

  \begin{equation} \label{kinetic_spont_gammad} \gamma_{d}=\frac{\frac{\sqrt{\pi}\left(3 V_{d}\right)^{1 / 3}}{S_{i}^{3 / 2}}\left(8 K C_{i} \sqrt{\pi S_{i}}-3 P \eta \right)-\sigma \alpha-\left(\sigma-\gamma_{m}\right) \beta}{\epsilon} \end{equation}
• CANCELLED solving for \(S_{i}\)
  
• DONE solving for \(V_{d}\)
  

  

  \begin{equation} \label{kinetic_spont_vd} V_{d}=\frac{\left(\frac{\left(\sigma \alpha+\left(\sigma-\gamma_{m}\right) \beta+\gamma_{d} \epsilon \right) S_{i}^{3 / 2}}{\left(8 K C_{i} \sqrt{\pi S_{i}}-3 P \eta \right) \sqrt{\pi}}\right)^{3}}{3} \end{equation}

configuration [12/13]

• import libraries, settings
  

  import napari
  import os
  from decimal import Decimal
  import math
  import numpy as np
  import pandas as pd
  import matplotlib.pyplot as plt  # plotting
  import matplotlib.animation as animation
  import mpl_toolkits.mplot3d.axes3d as axes3d
  import textwrap
  import seaborn as sns
  import scipy
  #sns.set(style='ticks')
  #sns.set_palette(sns.cubehelix_palette(10, start=.1, rot=-.75, reverse=True))
  sns.set_palette('hls',8)
  #sns.set_palette(sns.diverging_palette(250, 30, l=65, center="dark"))
  # sns.set_palette("coolwarm", 6)
  # from cycler import cycler
  SMALL_SIZE = 12
  MEDIUM_SIZE = 16
  BIGGER_SIZE = 20
  
  
  plt.rc('axes', titlesize=MEDIUM_SIZE)     # fontsize of the axes title
  plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
  plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
  plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
  plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
  plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title
  # default_cycler = (cycler(color=['indigo', 'darkgreen', 'purple', 'green', 'magenta', 'lime']))
  # plt.rc('axes', prop_cycle=default_cycler)
  

  #if kernel was started in notebook directory
  #fig_dir = 'figures/'
  #assuming kernel was started in home directory
  fig_dir = 'google_drive/grad_school/db_lab/code/analysis/llps_endocytosis_modeling/figures/'
  animation_dir = 'google_drive/grad_school/db_lab/code/analysis/llps_endocytosis_modeling/animation/'
  

  %matplotlib inline
  %gui qt
  import IPython
  from tabulate import tabulate
  
  class OrgFormatter(IPython.core.formatters.BaseFormatter):
      def __call__(self, obj):
          try:
              return tabulate(obj, headers='keys',
                              tablefmt='orgtbl', showindex='always')
          except:
              return None
  
  ip = get_ipython()
  ip.display_formatter.formatters['text/org'] = OrgFormatter()
  

• set default fixed parameter values
  
  • assigned here
    

    pi = np.pi
    vesicle_radius = 50. #nm
    vesicle_diameter = 2*vesicle_radius #nm
    drop_radius = 1500. #nm
    S_i = 4.*pi*vesicle_radius**2 #nm^2
    #H_i_range = np.linspace(0.001*vesicle_radius,2*vesicle_radius,1000) #nm
    h_range = np.linspace(0.001,1,1000) #dimensionless 0-1
    H_i_range = h_range*(S_i/pi)**(1/2)
    sigma = 150. #arbitrary, N/nm^2
    A_t = 1000 #nm, arbitrary
    K = 100. #arbitrary, N/nm
    C_i = 1/vesicle_radius #nm^-1
    gamma_d = 200. #arbitrary N/nm^2
    gamma_m = 175. #arbitrary, N/nm^2
    gamma_i = gamma_m
    V_d = 4.*pi*drop_radius**2 #nm^3
    P = 1000/(10**9) #N/nm^3
    
    params = {'vesicle_radius' : vesicle_radius,
        'vesicle_diameter' : vesicle_diameter,
        'drop_radius' : drop_radius,
        'S_i' : S_i,
        'h_range' : h_range,
        'H_i_range' : H_i_range,
        'sigma' : sigma,
        'A_t' : A_t,
        'K' : K,
        'C_i' : C_i,
        'gamma_d' : gamma_d,
        'gamma_m' : gamma_m,
        'gamma_i' : gamma_i,
        'V_d' : V_d,
        'P':P}
    

  • assigned from org table
    

    established in variables

    variable_df = pd.DataFrame(variable_list[1:], columns=variable_list[0])
    variable_df = variable_df.set_index('python variable')
    param_values = variable_df['value']
    variable_df
    
    

    vesicle_radius = param_values['vesicle_radius']
    vesicle_diameter = 2*vesicle_radius #nm
    drop_radius = param_values['drop_radius'] #nm
    S_i = 4.*pi*vesicle_radius**2 #nm^2
    #H_i_range = np.linspace(0.001*vesicle_radius,2*vesicle_radius,1000) #nm
    h_range = np.linspace(0.001,1,1000) #dimensionless 0-1
    H_i_range = h_range*(S_i/pi)**(1/2)
    sigma = param_values['sigma'] #arbitrary, N/nm^2
    A_t = param_values['A_t'] #nm, arbitrary
    K = param_values['K'] #arbitrary, N/nm
    C_i = 1/vesicle_radius #nm^-1
    gamma_d = param_values['gamma_d'] #arbitrary N/nm^2
    gamma_m = param_values['gamma_m'] #arbitrary, N/nm^2
    gamma_i = gamma_m
    #V_d = 4.*pi*drop_radius**2 #nm^3
    V_d = 4/3*pi*drop_radius**3 #nm^3
    P = param_values['P'] #N/nm^2
    
    params = {'vesicle_radius' : vesicle_radius,
        'vesicle_diameter' : vesicle_diameter,
        'drop_radius' : drop_radius,
        'S_i' : S_i,
        'h_range' : h_range,
        'H_i_range' : H_i_range,
        'sigma' : sigma,
        'A_t' : A_t,
        'K' : K,
        'C_i' : C_i,
        'gamma_d' : gamma_d,
        'gamma_m' : gamma_m,
        'gamma_i' : gamma_i,
        'V_d' : V_d,
        'P':P}
    

  • complemented from continuum model
    

    this is a bit adjusted so that there is a barrier due to the droplet

    N=10
    coat_radius = 50
    #coat_area = ((3*N)/(3*N+2))*4*np.pi*coat_radius**2
    coat_area = 4*np.pi*coat_radius**2
    total_radius = np.sqrt(coat_area/np.pi)*3
    Vd = 1/2*np.pi*200**3
    
    model_params = {
        #'time':40000,
        'time':1000000,
        'dt':1,
        'print_interval':10000,
        'write_interval':1000,
        'shape':'final',
        'P':0.0001,
        'sigma':0.01,
        'gamma_1':0,
        'gamma_2':0.05,
        'gamma_3':0.1,
        'gamma_d':0.2,
        'h1':0.00001,
        'h2':0.00001,
        'h3':0.00001,
        'H1':1/coat_radius,
        'K1':50,
        'H3':0,
        'K3':1,
        'polar_type':'integral',
        'G':1,
        'cell_wall':False,
        'N':N,
        'A_t':total_radius,
        'A_c':coat_area,
        'Vd':Vd,
        'emap':False,
        'tip_pull':0
    }
    
    V_d_range = np.linspace(5E6,5E7,100)
    

• define composite variables
  

  This is just to check what they are, I re-define them whenever I solve for F

  alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
  gamma_ratio = gamma_m/gamma_d
  alpha_refactored = np.cbrt((pi*(gamma_ratio+1)**3)/((gamma_ratio-1)*(gamma_ratio+2)**2))
  beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
  #epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
  epsilon = np.cbrt((8*pi)/((2+gamma_ratio)**2*(1-gamma_ratio)))
  eta = (S_i**(3/2))/(6*pi**(1/2))
  C = sigma*pi*A_t**2+S_i*(gamma_m-gamma_i)
  
  alpha, alpha_refactored, beta, epsilon, eta, C
  

• custom functions [12/13]
  
  • cube root
    

    def cube_root(i):
        return math.pow(abs(i),float(1)/3) * (1,-1)[i<0]
    

  • model functions
    
    • pre-analytical version
      
      • [ ] rewrite each function variable as dictionary key (outpus['theta_d'] instead of theta_d for the first one) and just return outputs at end
        • should eliminate some redundancy

      def llps_model(H_i=H_i_range, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                    C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d,
                    P=P):
          #H_i = h*(S_i/pi)**(1/2)
          # might need to be -gamma_m/gamma_d
          theta_d = np.arccos(gamma_m/gamma_d) #equation 8
          #to make a full array out of the constant theta_d value
          #theta_d = np.full_like(H_i, theta_d)
          A_i = ((S_i/pi)-H_i**2)**(1/2) #equation 7
          V_i = (1/6)*pi*H_i*(3*A_i**2+H_i**2) #equation 6
      # this is actually wrong!
      #    R_i = (A_i**2-H_i**2)/(2*H_i) #equation 10
          R_i = S_i/(2*pi*H_i) #equation 10 re-do
      # this is wrong!
      #    theta_i = np.arcsin(A_i/R_i) #equation 12
          theta_i = np.arccos(1-S_i/(2*pi*R_i**2)) #equation 12 re-do
          V_b = V_d+V_i # equation 5
          R_d_list = []
          #for V in V_b:
          #    R_d_list.append(cube_root((3*V)/(pi*(2+np.cos(theta_d))*(1-np.cos(theta_d))**2)))
          R_d = ((3*V_b)/(pi*(2+np.cos(theta_d))*(1-np.cos(theta_d))**2))**(1/3) #equation 4
          # R_d = np.array(R_d_list)
          S_d = 2*pi*R_d**2*(1-np.cos(theta_d)) #equation 11
          A_d = R_d*np.sin(theta_d) #equation 3
          S_m = pi*(A_d**2-A_i**2) #equation 9
          S_o = pi*(A_t**2-A_d**2) #equation 2
          membrane_tension = sigma*(S_o+S_m+S_i)
          bending_energy = (K/2)*S_i*((2/R_i)-(2*C_i))**2
          surface_tension = gamma_d*S_d
          wetting_energy = gamma_m*S_m+gamma_i*S_i
          droplet_energy = surface_tension - wetting_energy
          turgor_pressure = P*V_i
          F = membrane_tension+bending_energy+surface_tension-wetting_energy+turgor_pressure #equation 1
          non_droplet_energy = membrane_tension+bending_energy+turgor_pressure
      
          outputs = {'theta_d':theta_d,'A_i':A_i,'V_i':V_i, 'R_i':R_i,
                     '2÷R_i':2/R_i,'theta_i':theta_i,'V_b':V_b,'R_d':R_d,
                    'S_d':S_d,'A_d':A_d,'S_m':S_m,'S_o':S_o,
                     'F':F,'membrane_tension':membrane_tension,
                     'bending_energy':bending_energy, 'bending_energy_flipped':-bending_energy,
                     'surface_tension':surface_tension,
                     'wetting_energy':-wetting_energy, 'wetting_energy_flipped':wetting_energy,
                     'droplet_energy':droplet_energy, 'turgor_pressure':turgor_pressure,
                     'non_droplet_energy':non_droplet_energy}
          return outputs
      

      plt.plot(H_i_range/100,llps_model()['F'])
      plt.xlabel('h')
      plt.ylabel('F')
      plt.tight_layout()
      plt.savefig(fig_dir+'llps_model_test.png')
      

    • post-analytical version
      

      def llps_model_recbrt(h=h_range, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                     C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i,
                     V_d=V_d, P=P):
          
          alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
          eta = (S_i**(3/2))/(6*pi**(1/2))
          C = sigma*pi*A_t**2+S_i*(gamma_m-gamma_i)
          
          F = ((sigma*alpha+(sigma-gamma_m)*beta+gamma_d*epsilon)* np.cbrt(((S_i**(3/2))/(pi**(1/2)))*(3*h/2-h**3)+3*V_d)**2+ (sigma-gamma_m)*S_i*h**2+ 2*K*(C_i*S_i**(1/2)-2*pi**(1/2)*h)**2+ P*eta*(3*h-2*h**3)+ C )
      
      #    print(alpha, beta, epsilon, eta, C)
      
          return F
      

      plt.plot(h_range,llps_model_recbrt())
      plt.xlabel('h')
      plt.ylabel('F')
      plt.tight_layout()
      plt.savefig(fig_dir+'llps_model_recbrt_test.png')
      

  • general plotting function
    
    • figure out some way to get the right units for each thing being plotted

    def plot_simulation(x=H_i_range/vesicle_diameter, output='F',
                     H_i=H_i_range, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                     C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P,
                     label=None, variable_df=variable_df):
        outputs = llps_model(H_i=H_i_range, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                  C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
        y = outputs[output]
        plt.plot(x,y,label=label)
        plt.ylabel(variable_df.loc[output]['description']+' ('+variable_df.loc[output]['units']+')')
    

    • testing
      

      plot_simulation()
      

  • normalized plots
    

    def plot_simulation_normalized(x=H_i_range/vesicle_diameter, output='F',
                    H_i=H_i_range, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                    C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P,
                                    label=None, color=None, linestyle='-', linewidth=None):
        outputs = llps_model(H_i=H_i_range, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d)
        y_sub = outputs[output]-np.min(outputs[output])
        y_norm = y_sub/np.max(y_sub)
        plt.plot(x,y_norm,label=label, color=color, linestyle=linestyle, linewidth=linewidth)
        plt.ylabel(output)
        plt.xlabel('h: invagination distance/vesicle diameter')
    

  • delta plots
    

    def plot_simulation_delta(x=H_i_range/vesicle_diameter, output='F',
                     H_i=H_i_range, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                     C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P,
                     label=None, variable_df=variable_df, color=None, linestyle='-', linewidth=None):
        outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                  C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
        output_0 = llps_model(H_i=H_i[0], sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                  C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
        y = outputs[output]-output_0[output]
        plt.plot(x,y,label=label, color=color,linestyle=linestyle, linewidth=linewidth)
        #plt.ylabel('∆ '+variable_df.loc[output]['description']+' ('+variable_df.loc[output]['units']+')')
        plt.xlabel('h: invagination distance/vesicle diameter')
    

  • show geometry for given invagination distance
    

    def show_geometry(H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                      C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P,
                      fig_width=10, fig_height=5,
                      num_plots=1, top=20, bottom=-250, left=-500, right=500,
                      intermediate=False):
        outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                 C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
    
        A_d = outputs['A_d']
        A_i = outputs['A_i']
        theta_d = outputs['theta_d']
        theta_i = outputs['theta_i']
        R_d = outputs['R_d']
        R_i = outputs['R_i']
    
        #fig, ax = plt.subplots(1,num_plots)
        #plt.subplot(1, num_plots, 1)
    
        fig = plt.figure()
        fig.set_size_inches((fig_width,fig_height))
        ax = fig.add_subplot(1, num_plots, 1)
    
        plt.hlines([0,0],[-A_t,A_d],[-A_d,A_t], color='r', label='$S_o$')
        plt.hlines([0,0],[-A_d,A_i],[-A_i,A_d], color='g', label='$S_m$')
    
        if V_d > 0:
          points_d = np.linspace(pi/2-theta_d, theta_d+pi/2, 1000)
          x_d = R_d*np.cos(points_d)
          y_d = R_d*np.sin(points_d)
          y_d_adj = y_d-np.min(y_d)
          plt.plot(x_d,-y_d_adj, color='b', label='$S_d$')
    
        points_i = np.linspace(pi/2-theta_i, theta_i+pi/2, 1000)
        x_i = R_i*np.cos(points_i)
        y_i = R_i*np.sin(points_i)
        y_i_adj = y_i-np.min(y_i)
        plt.plot(x_i,-y_i_adj, color='purple', label='$S_i$')
    
        ax.set_ylim(top=top, bottom=bottom)
        ax.set_xlim(left=left, right=right)
        ax.set_aspect('equal')
        ax.set(xlabel = 'distance along membrane (nm)',
               ylabel = 'distance into cytoplasm (nm)')
        plt.legend()
        plt.tight_layout()
    
        if intermediate==True:
               return fig, ax
    

    • testing
      

      show_geometry(H_i=100)
      

  • get geometry to feed into continuum model
    

    def get_geometry(H_i, N, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                      C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P):
        outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                 C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
    
        r = np.linspace(0, A_t, 3*N+1)
        z = np.full(3*N+1,0)
    
        A_d = outputs['A_d']
        A_i = outputs['A_i']
        theta_d = outputs['theta_d']
        theta_i = outputs['theta_i']
        R_d = outputs['R_d']
        R_i = outputs['R_i']
    
        r[N:2*N+1] = np.linspace(A_i,A_d,N+1)
        r[2*N:3*N+1] = np.linspace(A_d,A_t,N+1)
    
        points_i = np.linspace(pi/2, pi/2-theta_i, N+1)
        x_i = R_i*np.cos(points_i)
        y_i = R_i*np.sin(points_i)
        y_i_adj = y_i-np.min(y_i)
        r[0:N+1] = x_i
        z[0:N+1] = -y_i_adj
    
        return(r,z)
    

    • testing
      

      show_geometry(H_i=100)
      

  • 3D geometry
    
    • inspired by noauthor_python_nodate
    • need to turn this into a function
    • points
      

      %gui qt
      viewer = napari.Viewer(ndisplay=3)
      
      spherecoords = []
      for H_i in np.linspace(1,100,100):
          outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K, C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
      
          A_d = outputs['A_d']
          A_i = outputs['A_i']
          theta_d = outputs['theta_d']
          theta_i = outputs['theta_i']
          R_d = outputs['R_d']
          R_i = outputs['R_i']
          theta, phi = np.linspace(0,2* np.pi, 200), np.linspace(pi-theta_d, np.pi, 200)
          THETA, PHI = np.meshgrid(theta, phi)
          R = R_d
          X = R * np.sin(PHI) * np.cos(THETA)
          Y = R * np.sin(PHI) * np.sin(THETA)
          Z = R * np.cos(PHI)
          Z = Z-np.max(Z)
          height = np.full_like(Z, H_i)
          spherecoords.append(np.array([np.ravel(height), np.ravel(Z), np.ravel(Y), np.ravel(X)]).T)
      
      spherecoords = np.concatenate(spherecoords)
      viewer.add_points(np.array(spherecoords), name='droplet', size=20, face_color='cyan',
                        edge_color = 'transparent', blending='additive', opacity=0.6)
      
      spherecoords = []
      for H_i in np.linspace(1,100,100):
          outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K, C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
      
          A_d = outputs['A_d']
          A_i = outputs['A_i']
          theta_d = outputs['theta_d']
          theta_i = outputs['theta_i']
          R_d = outputs['R_d']
          R_i = outputs['R_i']
          theta, phi = np.linspace(0,2* np.pi, 200), np.linspace(pi-theta_i, np.pi, 200)
          THETA, PHI = np.meshgrid(theta, phi)
          R = R_i
          X = R * np.sin(PHI) * np.cos(THETA)
          Y = R * np.sin(PHI) * np.sin(THETA)
          Z = R * np.cos(PHI)
          Z = Z-np.max(Z)
          height = np.full_like(Z, H_i)
          spherecoords.append(np.array([np.ravel(height), np.ravel(Z), np.ravel(Y), np.ravel(X)]).T)
      
      spherecoords = np.concatenate(spherecoords)
      viewer.add_points(np.array(spherecoords), name='coat', size=20, face_color='purple',
                        edge_color = 'transparent', blending='additive', opacity=0.6)
      
      spherecoords = []
      for H_i in np.linspace(1,100,100):
          outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K, C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
      
          A_d = outputs['A_d']
          A_i = outputs['A_i']
          theta_d = outputs['theta_d']
          theta_i = outputs['theta_i']
          theta_m = 0.00000001
          R_d = outputs['R_d']
          R_i = outputs['R_i']
          theta, phi = np.linspace(0,2* np.pi, 200), np.linspace(pi-theta_m, np.pi, 200)
          THETA, PHI = np.meshgrid(theta, phi)
          R = 0.1*A_t/np.sin(theta_m)
          X = R * np.sin(PHI) * np.cos(THETA)
          Y = R * np.sin(PHI) * np.sin(THETA)
          Z = R * np.cos(PHI)
          Z = Z-np.max(Z)+0
          height = np.full_like(Z, H_i)
          spherecoords.append(np.array([np.ravel(height), np.ravel(Z), np.ravel(Y), np.ravel(X)]).T)
      
      spherecoords = np.concatenate(spherecoords)
      viewer.add_points(np.array(spherecoords), name='cell wall', size=20, face_color='red',
                        edge_color = 'transparent', blending='additive', opacity=0.6)
      
      

      • I don’t understand why this doesn’t really work

      vertices = spherecoords
      #faces = np.array([[0, 0, 1, 2], [0, 1, 2, 3], [0, 2, 3, 4]])
      faces_list = []
      counter = 0
      for H_i in range(len(H_i_range)):
          for point in range(38):
              pointa = counter + point
              faces_list.append([pointa, pointa+1, pointa+2])
          counter += 40
      faces = np.array(faces_list)
      values = np.linspace(0, 1, len(vertices))
      surface = (vertices, faces, values)
      
      viewer.add_surface(surface, colormap='viridis')  # add the surface
      

    • surface
      

      examples:

      vertices = np.array([[0, 0, 0, 0], [0, 10, 0, 20], [0, 20, 10, 0], [0, 0, 10, 10],
                           [1, 20, 0, 0], [1, 10, 0, 20], [1, 20, 10, 0], [1, 0, 10, 10]])
      faces = np.array([[0, 1, 2], [1, 2, 3],
                        [4,5,6],[5,6,7]])
      values = np.linspace(0.9, 1, len(vertices))
      surface = (vertices, faces, values)
      
      viewer = napari.view_surface(surface)  # add the surface
      
      facesloop = []
      
      for t in (0,1):
          vert = vertices[vertices[:,0] == t]
          facest = scipy.spatial.Delaunay(vert[:,2:]).simplices + len(vert)*t
          facesloop.append(facest)
      
      facesd = np.concatenate(facesloop)
      
      viewer = napari.view_surface((vertices, facesd, values))
      

      coords = np.array([np.ravel(height), np.ravel(Z), np.ravel(Y), np.ravel(X)]).T[:,1:]
      faces = scipy.spatial.Delaunay( coords[:,1:]).simplices
      values = np.linspace(0,1,len(coords))
      viewer = napari.view_surface((coords[:,1:], faces, values))
      

      %gui qt
      viewer = napari.Viewer(ndisplay=3)
      
      spherecoords = []
      faces = []
      for H_i in range(1,101):
          outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K, C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
      
          A_d = outputs['A_d']
          A_i = outputs['A_i']
          theta_d = outputs['theta_d']
          theta_i = outputs['theta_i']
          R_d = outputs['R_d']
          R_i = outputs['R_i']
          theta, phi = np.linspace(0,2* np.pi, 50), np.linspace(pi-theta_d, np.pi, 50)
          THETA, PHI = np.meshgrid(theta, phi)
          R = R_d
          X = R * np.sin(PHI) * np.cos(THETA)
          Y = R * np.sin(PHI) * np.sin(THETA)
          Z = R * np.cos(PHI)
          Z = Z-np.max(Z)
          height = np.full_like(Z, H_i)
          vertices = np.array([np.ravel(height), np.ravel(Z), np.ravel(Y), np.ravel(X)]).T
          spherecoords.append(vertices)
      
          #face = scipy.spatial.Delaunay(vertices[:,2:]).simplices + len(vertices)*(H_i-1)
          if H_i == 1:
              face = scipy.spatial.Delaunay(vertices[:,2:]).simplices + len(vertices)*(H_i-1)
          else:
              face = faces[-1] + len(vertices)
          faces.append(face)
      
      spherecoords = np.concatenate(spherecoords)
      faces = np.concatenate(faces)
      surface_values = np.full(len(spherecoords),1)
      viewer.add_surface((spherecoords,faces,surface_values), name='droplet', colormap='green',
                         blending='translucent', shading='smooth')
      
      spherecoords = []
      faces = []
      for H_i in range(1,101):
          outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K, C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
      
          A_d = outputs['A_d']
          A_i = outputs['A_i']
          theta_d = outputs['theta_d']
          theta_i = outputs['theta_i']
          R_d = outputs['R_d']
          R_i = outputs['R_i']
          theta, phi = np.linspace(0,2* np.pi, 50), np.linspace(pi-theta_i, np.pi, 50)
          THETA, PHI = np.meshgrid(theta, phi)
          R = R_i
          X = R * np.sin(PHI) * np.cos(THETA)
          Y = R * np.sin(PHI) * np.sin(THETA)
          Z = R * np.cos(PHI)
          Z = Z-np.max(Z)
          height = np.full_like(Z, H_i)
          vertices = np.array([np.ravel(height), np.ravel(Z), np.ravel(Y), np.ravel(X)]).T
          spherecoords.append(vertices)
      
          if H_i == 1:
              face = scipy.spatial.Delaunay(vertices[:,2:]).simplices + len(vertices)*(H_i-1)
          else:
              face = faces[-1] + len(vertices)
          faces.append(face)
      
      spherecoords = np.concatenate(spherecoords)
      faces = np.concatenate(faces)
      surface_values = np.full(len(spherecoords),1)
      viewer.add_surface((spherecoords,faces,surface_values), name='coat', colormap='bop blue',
                         blending='additive', shading='smooth')
      
      spherecoords = []
      faces = []
      for H_i in range(1,101):
          outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K, C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
      
          A_d = outputs['A_d']
          A_i = outputs['A_i']
          theta_d = outputs['theta_d']
          theta_i = outputs['theta_i']
          theta_m = 0.00000001
          R_d = outputs['R_d']
          R_i = outputs['R_i']
          theta, phi = np.linspace(0,2* np.pi, 50), np.linspace(pi-theta_m, np.pi, 50)
          THETA, PHI = np.meshgrid(theta, phi)
          R = 0.1*A_t/np.sin(theta_m)
          X = R * np.sin(PHI) * np.cos(THETA)
          Y = R * np.sin(PHI) * np.sin(THETA)
          Z = R * np.cos(PHI)
          Z = Z-np.max(Z)+0
          height = np.full_like(Z, H_i)
          vertices = np.array([np.ravel(height), np.ravel(Z), np.ravel(Y), np.ravel(X)]).T
          spherecoords.append(vertices)
      
          #face = scipy.spatial.Delaunay(vertices[:,2:]).simplices + len(vertices)*(H_i-1)
          if H_i == 1:
              face = scipy.spatial.Delaunay(vertices[:,2:]).simplices + len(vertices)*(H_i-1)
          else:
              face = faces[-1] + len(vertices)
          faces.append(face)
      
      spherecoords = np.concatenate(spherecoords)
      faces = np.concatenate(faces)
      surface_values = np.full(len(spherecoords),1)
      viewer.add_surface((spherecoords,faces,surface_values), name='cell wall', colormap='bop blue',
                         blending='opaque', shading='smooth')
      
      def update_frame():
          """Update frame or time"""
          viewer.text_overlay.text = 'h = %1.2f\nblue: membrane\ngreen: condensate' % ((viewer.dims.current_step[0]+1)/100)
      
      viewer.dims.events.current_step.connect(lambda event: update_frame())
      
      viewer.text_overlay.visible=True
      viewer.text_overlay.font_size=20
      viewer.text_overlay.color = 'white'
      viewer.scale_bar.unit = 'nm'
      viewer.scale_bar.visible=True
      viewer.scale_bar.font_size=10
      

      animation script

      from napari_animation import Animation
      
      animation = Animation(viewer)
      viewer.reset_view()
      viewer.dims.current_step = (0, 0,0,0)
      animation.capture_keyframe(steps=15)
      animation.capture_keyframe(steps=15)
      viewer.camera.zoom = 2.6
      viewer.camera.angles=(-0.026467524771271397, 0.5719902891444765, 11.550928988495166)
      animation.capture_keyframe(steps=50)
      viewer.dims.current_step = (viewer.dims.range[0][1]-1, 0,0,0)
      animation.capture_keyframe(steps=int(viewer.dims.range[0][1]-1))
      viewer.reset_view()
      animation.capture_keyframe(steps=50)
      viewer.dims.current_step = (0, 0,0,0)
      animation.capture_keyframe(steps=int(viewer.dims.range[0][1]-1))
      animation.animate('~/google_drive/grad_school/db_lab/code/analysis/llps_endocytosis_modeling/scripted_surface_animation.mov')
      

    • probably don’t need this stuff
      

      num_plots = 1
      outputs = llps_model(H_i=90, sigma=sigma, A_t=A_t, S_i=S_i, K=K, C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
      
      A_d = outputs['A_d']
      A_i = outputs['A_i']
      theta_d = outputs['theta_d']
      theta_i = outputs['theta_i']
      R_d = outputs['R_d']
      R_i = outputs['R_i']
      
      #fig, ax = plt.subplots(1,num_plots)
      #plt.subplot(1, num_plots, 1)
      
      #fig = plt.figure()
      #fig.set_size_inches((fig_width,fig_height))
      #ax = fig.add_subplot(1, num_plots, 1)
      
      #plt.hlines([0,0],[-A_t,A_d],[-A_d,A_t], color='r', label='$S_o$')
      #plt.hlines([0,0],[-A_d,A_i],[-A_i,A_d], color='g', label='$S_m$')
      
      if V_d > 0:
          points_d = np.linspace(pi/2-theta_d, theta_d+pi/2, 1000)
          x_d = R_d*np.cos(points_d)
          y_d = R_d*np.sin(points_d)
          y_d_adj = y_d-np.min(y_d)
      
      points_i = np.linspace(pi/2-theta_i, theta_i+pi/2, 1000)
      x_i = R_i*np.cos(points_i)
      y_i = R_i*np.sin(points_i)
      y_i_adj = y_i-np.min(y_i)
      
      xy_i = np.array([y_i,x_i]).T
      viewer = napari.Viewer()
      viewer.add_points(xy_i)
      

      nx, ny, nz = (3, 3, 4)
      x = np.linspace(0, 10, nx)
      y = np.linspace(0, 10, ny)
      z = np.linspace(0, 10, nz)
      
      xv, yv, zv = np.meshgrid(x, y, z)
      zyxv = np.array([np.ravel(zv), np.ravel(yv), np.ravel(xv)]).T
      zyxv
      
      viewer = napari.Viewer()
      viewer.add_points(zyxv)
      

      import numpy as np
      import matplotlib.pyplot as plt
      from mpl_toolkits import mplot3d
      
      ax = plt.axes(projection="3d")
      xlist=np.linspace(-1.0,1.0,10)*100
      ylist=np.linspace(-1.0,1.0,10)*100
      r=np.linspace(1.0,1.0,10)*100
      X,Y= np.meshgrid(xlist,ylist)
      
      Z=np.sqrt(r**2-X**2-Y**2) #Use np.sqrt like you had before
      
      spherecoords = np.array([np.ravel(Z), np.ravel(Y), np.ravel(X)]).T
      viewer.add_points(spherecoords)
      
      vertices = spherecoords
      faces = np.array([[0, 50, 99], [1, 2, 3]])
      values = np.linspace(1, 1, len(vertices))
      surface = (vertices, faces, values)
      
      viewer.add_surface(surface)
      
      cp=ax.plot_wireframe(X,Y,Z,color="r")
      cp=ax.plot_wireframe(X,Y,-Z,color="r") # Now plot the bottom half
      
      
      plt.title('2D Contour Plot of The unit sphere')
      plt.show()
      

  • side-by-side delta plot and equilibrium geometry
    

    def plot_simulation_geometry(x=H_i_range/vesicle_diameter, output='F',
                     H_i=H_i_range, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                     C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P,
                     fig_width=15, fig_height=5, label=None, variable_df=variable_df,
                    num_plots=2, top=20, bottom=-250, left=-500, right=500):
        outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                  C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
        min_energy_ix = np.argmin(outputs['F'])
        H_i_eq = H_i[min_energy_ix]
        fig, ax1 = show_geometry(H_i=H_i_eq, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                       C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P,
                       fig_width=fig_width, fig_height=fig_height,
                       num_plots=num_plots, top=top, bottom=bottom, left=left, right=right,
                       intermediate=True)
        plt.title('geometry at equilibrium')
        ax2 = fig.add_subplot(1,2,2)
        #y = outputs[output]
        #plt.plot(x,y,label=label)
        #plt.ylabel(output)
    
        plot_simulation_delta(x=x, output=output,
                     H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                     C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P,
                     label=label, variable_df=variable_df)
    
        plt.tight_layout()
    

    • testing
      

      plot_simulation_geometry()
      

  • animated side-by-side geometry and energy plot, scan volume
    

    V_d_range=np.linspace(V_d*0.1,V_d*10,100)
    def animate_simulation_geometry(x=H_i_range/vesicle_diameter, output='F',
                    H_i=H_i_range, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                    C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d_range=V_d_range, P=P,
                    fig_width=15, fig_height=5, label=None, color=None, variable_df=variable_df,
                    geometry_top=20, geometry_bottom=-250, geometry_left=-500, geometry_right=500,
                    plot_top=20, plot_bottom=-250, plot_left=-500, plot_right=500):
        i = 0
        V_d=V_d_range[i]
    
        fig, (ax_geometry, ax_plot) = plt.subplots(nrows=1, ncols=2, figsize=(15,5))
    
        outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                    C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
        output_0 = llps_model(H_i=H_i[0], sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
        min_energy_ix = np.argmin(outputs['F'])
        H_i_eq = H_i[min_energy_ix]
        outputs_eq = llps_model(H_i=H_i_eq, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
    
        A_d = outputs_eq['A_d']
        A_i = outputs_eq['A_i']
        theta_d = outputs_eq['theta_d']
        theta_i = outputs_eq['theta_i']
        R_d = outputs_eq['R_d']
        R_i = outputs_eq['R_i']
    
        top = geometry_top
        bottom = geometry_bottom
        left = geometry_left
        right = geometry_right
    
        S_o_line_l, = ax_geometry.plot([-A_t,-A_d],[0,0], color='b', label='$S_o$')
        S_o_line_r, = ax_geometry.plot([A_d,A_t],[0,0], color='b')
        S_m_line_l, = ax_geometry.plot([-A_d,-A_i],[0,0], color='cyan', label='$S_m$')
        S_m_line_r, = ax_geometry.plot([A_i,A_d],[0,0], color='cyan')
        if V_d > 0:
            points_d = np.linspace(pi/2-theta_d, theta_d+pi/2, 1000)
            x_d = R_d*np.cos(points_d)
            y_d = R_d*np.sin(points_d)
            y_d_adj = y_d-np.min(y_d)
            S_d_line, = ax_geometry.plot(x_d,-y_d_adj, color='g', label='$S_d$')
            if geometry_bottom == None:
                bottom = y_d_adj.max()*(-1.1)
            if geometry_left == None:
                left = x_d.max()*(-1.2)
            if geometry_right == None:
                right = x_d.max()*(1.2)
        points_i = np.linspace(pi/2-theta_i, theta_i+pi/2, 1000)
        x_i = R_i*np.cos(points_i)
        y_i = R_i*np.sin(points_i)
        y_i_adj = y_i-np.min(y_i)
        S_i_line, = ax_geometry.plot(x_i,-y_i_adj, color='purple', label='$S_i$')
    
        ax_geometry.set_ylim(top=top, bottom=bottom)
        ax_geometry.set_xlim(left=left, right=right)
        ax_geometry.set_aspect('equal')
        ax_geometry.set(xlabel = 'distance along membrane (nm)',
                ylabel = 'distance into cytoplasm (nm)')
        ax_geometry.legend()
        ax_geometry.set_title('geometry at equilibrium')
    
    
        surface_tension = outputs['surface_tension']-output_0['surface_tension']
        wetting_energy = outputs['wetting_energy']-output_0['wetting_energy']
        fd = outputs['droplet_energy']-output_0['droplet_energy']
        plot_st, = ax_plot.plot(x,surface_tension,label='$\gamma_d S_d$', color='green')
        plot_we, = ax_plot.plot(x,wetting_energy,label='$\gamma_m S_m$', color='cyan')
        plot_fd, = ax_plot.plot(x,fd,label='$F_d$', color='black')
        plt.hlines(0,0,1,linestyles=':',colors='gray')
        if plot_top == None:
            ax_plot.set_ylim(top=surface_tension.max()*1.2,
                            bottom=wetting_energy.min()-surface_tension.max()*0.2)
        else:
            ax_plot.set_ylim(top=plot_top, bottom=plot_bottom)
        ax_plot.legend()
        ax_plot.set(xlabel='h: invagination distance/vesicle diameter',
            ylabel='∆ energy $(pN \cdot nm)$')
    
        plt.tight_layout()
    
        def animate(i):
            V_d=V_d_range[i]
    
            outputs = llps_model(H_i=H_i, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                        C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
            output_0 = llps_model(H_i=H_i[0], sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                    C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
            min_energy_ix = np.argmin(outputs['F'])
            H_i_eq = H_i[min_energy_ix]
            outputs_eq = llps_model(H_i=H_i_eq, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                    C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i, V_d=V_d, P=P)
    
            A_d = outputs_eq['A_d']
            A_i = outputs_eq['A_i']
            theta_d = outputs_eq['theta_d']
            theta_i = outputs_eq['theta_i']
            R_d = outputs_eq['R_d']
            R_i = outputs_eq['R_i']
    
            top = geometry_top
            bottom = geometry_bottom
            left = geometry_left
            right = geometry_right
    
            #S_o_line = ax_geometry.hlines([0,0],[-A_t,A_d],[-A_d,A_t], color='r', label='$S_o$')
            S_o_line_l.set_data([-A_t,-A_d],[0,0])
            S_o_line_r.set_data([A_d,A_t],[0,0])
            #S_m_line = ax_geometry.hlines([0,0],[-A_d,A_i],[-A_i,A_d], color='g', label='$S_m$')
            S_m_line_l.set_data([-A_d,-A_i],[0,0])
            S_m_line_r.set_data([A_i,A_d],[0,0])
            if V_d > 0:
                points_d = np.linspace(pi/2-theta_d, theta_d+pi/2, 1000)
                x_d = R_d*np.cos(points_d)
                y_d = R_d*np.sin(points_d)
                y_d_adj = y_d-np.min(y_d)
                S_d_line.set_data(x_d,-y_d_adj)
                if geometry_bottom == None:
                    bottom = y_d_adj.max()*(-1.1)
                if geometry_left == None:
                    left = x_d.max()*(-1.2)
                if geometry_bottom == None:
                    right = x_d.max()*(1.2)
            points_i = np.linspace(pi/2-theta_i, theta_i+pi/2, 1000)
            x_i = R_i*np.cos(points_i)
            y_i = R_i*np.sin(points_i)
            y_i_adj = y_i-np.min(y_i)
            S_i_line.set_data(x_i,-y_i_adj)
            ax_geometry.set_ylim(top=top, bottom=bottom)
            ax_geometry.set_xlim(left=left, right=right)
            ax_geometry.set_aspect('equal')
    
            surface_tension = outputs['surface_tension']-output_0['surface_tension']
            wetting_energy = outputs['wetting_energy']-output_0['wetting_energy']
            fd = outputs['droplet_energy']-output_0['droplet_energy']
            plot_st.set_data(x,surface_tension)
            plot_we.set_data(x,wetting_energy)
            plot_fd.set_data(x,fd)
            if plot_top == None:
                ax_plot.set_ylim(top=surface_tension.max()*1.2,
                                bottom=wetting_energy.min()-surface_tension.max()*0.2)
    
            plt.tight_layout()
        ani = animation.FuncAnimation(
            fig, animate, frames=len(V_d_range), interval=100)
        ani.save(os.path.join(animation_dir, 'vdscan.mp4'), dpi=300)
    
        plt.show()
    

  • side-by-side F plots and parameter phase diagram
    

    def boundary_test(x_param, y_param, x_range, boundary, params=params,
                      factors=[0.7,1,1.3], y_offset=0.2, x_scale='log'):
        fparams = params.copy()
        fig, (ax1, ax2)  = plt.subplots(nrows=1, ncols=2, figsize=(10,4))
    
        if x_scale == 'log':
            ax1.semilogx(x_range, boundary, color='black')
        else:
            ax1.plot(x_range, boundary, color='black')
    
        for factor in factors:
            middle = len(x_range)//2
            x_pick = x_range[middle]*factor
            boundary_pick = boundary[middle]+(factor-1)*np.sqrt(boundary[middle]**2)*y_offset
            fparams[x_param] = x_pick
            fparams[y_param] = boundary_pick
            F = llps_model_recbrt(h=fparams['h_range'], sigma=fparams['sigma'], A_t=fparams['A_t'],
                    S_i=fparams['S_i'], K=fparams['K'], C_i=fparams['C_i'], gamma_d=fparams['gamma_d'],
                    gamma_m=fparams['gamma_m'], gamma_i=fparams['gamma_i'], V_d=fparams['V_d'], P=fparams['P'])
            delta_F = F-F[0]
    
            ax1.scatter(x_pick, boundary_pick)
            ax2.plot(h_range,delta_F)
    
        ax2.hlines(0, 0,1.05, color='black')
    
        ax1.set(xlabel=x_param, ylabel=y_param)
        ax2.set(xlabel='h', ylabel='∆F')
        plt.tight_layout()
    

  • global stability parameter analysis [4/4]
    
    • DONE \(K\)
      

      from therm_spont_k

      def thermo_spont_K(sigma=sigma, gamma_m=gamma_m, gamma_d=gamma_d,
                        V_d=V_d, S_i=S_i, P=P, C_i=C_i):
          
          alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
          eta = (S_i**(3/2))/(6*pi**(1/2))
          
          K = ((-(sigma*alpha+(sigma-gamma_m)*beta+gamma_d*epsilon)*
                (np.cbrt(3*V_d+(S_i**(3/2))/(2*pi**(1/2)))**2-
                  np.cbrt(3*V_d)**2)+2*P*eta-(sigma-gamma_m)*S_i-
                3*P*eta)/(8*pi-8*C_i*(pi*S_i)**(1/2)))
          
          return K
      

      • testing
        

        V_d_range = np.logspace(6,9)
        K_boundary = thermo_spont_K(V_d=V_d_range)
        boundary_test(x_param='V_d', y_param='K', x_range=V_d_range, boundary=K_boundary)
        

        C_i_range = np.logspace(-4,-2.1,1000)
        K_boundary = thermo_spont_K(C_i=C_i_range)
        boundary_test(x_param='C_i', y_param='K', x_range=C_i_range, boundary=K_boundary,
                      y_offset=-0.3)
        

        • this kind of makes sense - the lower the spontaneous curvature (flatter), the more the energetics are sensitive to rigidity

        x_range = np.logspace(6,9)
        thermo_boundary = thermo_spont_K(V_d=x_range)
        kinetic_boundary = kinetic_spont_K(V_d=x_range)
        both_boundary_test(x_param='V_d', y_param='K', x_range=x_range,
                           thermo_boundary=thermo_boundary, kinetic_boundary=kinetic_boundary,
                           y_offset=2)
        

    • DONE \(P\)
      
      • from therm_spont_p

      def thermo_spont_P(sigma=sigma, gamma_m=gamma_m, gamma_d=gamma_d,
                        V_d=V_d, S_i=S_i, K=K, C_i=C_i):
          
          alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
          eta = (S_i**(3/2))/(6*pi**(1/2))
          
          P = ((sigma*alpha + (sigma - gamma_m)*beta + gamma_d*epsilon)*(np.cbrt(3*V_d)**2 -
               np.cbrt(S_i**(3/2)/(2*np.sqrt(pi)) + 3*V_d)**2) -
               (sigma-gamma_m)*S_i - 8*K*(pi - C_i*np.sqrt(pi*S_i)))/ eta
          
          return P
      

      • testing
        

        x_range = np.logspace(6,9)
        boundary = thermo_spont_P(V_d=x_range)
        boundary_test(x_param='V_d', y_param='P', x_range=x_range, boundary=boundary,
                      y_offset=-0.3)
        

    • DONE \(\sigma\)
      
      • from therm_spont_sigma

      def thermo_spont_sigma(P=P, gamma_m=gamma_m, gamma_d=gamma_d,
                        V_d=V_d, S_i=S_i, K=K, C_i=C_i):
          
          alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
          eta = (S_i**(3/2))/(6*pi**(1/2))
      
          sigma = ((-gamma_m*S_i+8*K*(pi-C_i*np.sqrt(pi*S_i))+P*eta)/
                  ((3*V_d)**(2/3)-((S_i**(3/2))/
                      (2*np.sqrt(pi))+3*V_d)**(2/3))-
                   gamma_d*epsilon+gamma_m*beta)/(alpha+beta-
                  (S_i)/((3*V_d)**(2/3)-((S_i**(3/2))/
                          (2*np.sqrt(pi))+3*V_d)**(2/3)))
          
          return sigma
      

      • testing
        

        x_range = np.logspace(6,9)
        boundary = thermo_spont_sigma(V_d=x_range)
        boundary_test(x_param='V_d', y_param='sigma', x_range=x_range, boundary=boundary,
                      y_offset=-0.02)
        

    • DONE \(C_{i}\)
      
      • from therm_spont_ci

      def thermo_spont_C_i(P=P, gamma_m=gamma_m, gamma_d=gamma_d,
                        V_d=V_d, S_i=S_i, K=K, sigma=sigma):
          
          alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
          eta = (S_i**(3/2))/(6*pi**(1/2))
      
          C_i = (pi-((sigma*alpha+(sigma-gamma_m)*beta+gamma_d*epsilon)*(np.cbrt(3*V_d)**2-np.cbrt((S_i**(3/2))/(2*np.sqrt(pi))+3*V_d)**2)-(sigma-gamma_m)*S_i-P*eta)/(8*K))/(np.sqrt(pi*S_i))
          
          return C_i
      

      • testing
        

        x_range = np.logspace(2,7)
        boundary = thermo_spont_C_i(V_d=x_range)
        boundary_test(x_param='V_d', y_param='C_i', x_range=x_range, boundary=boundary,
                      y_offset=0.2)
        

        • it doesn’t make sense that C_i would ever be negative…
  • local stability parameter analysis [8/8]
    
    • DONE \(K\)
      

      from kinetic_spont_k

      def kinetic_spont_K(sigma=sigma, gamma_m=gamma_m, gamma_d=gamma_d,
                        V_d=V_d, S_i=S_i, P=P, C_i=C_i):
          
          alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
          eta = (S_i**(3/2))/(6*pi**(1/2))
          
          K = (((sigma*alpha+(sigma-gamma_m)*beta+gamma_d*epsilon)*S_i**(3/2))/(np.sqrt(pi)*np.cbrt(3*V_d)) + 3*P*eta)/(8*C_i*np.sqrt(pi*S_i))
          
          return K
      

      • testing
        

        x_range = np.logspace(6,9)
        boundary = kinetic_spont_K(V_d=x_range)
        boundary_test(x_param='V_d', y_param='K', x_range=x_range, boundary=boundary,
                      y_offset=0.5)
        

    • DONE \(P\)
      

      from kinetic_spont_p

      def kinetic_spont_P(sigma=sigma, gamma_m=gamma_m, gamma_d=gamma_d,
                        V_d=V_d, S_i=S_i, K=K, C_i=C_i):
          
          alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
          eta = (S_i**(3/2))/(6*pi**(1/2))
          
          P = (-((sigma*alpha + (sigma - gamma_m)*beta +gamma_d*epsilon)*S_i**(3/2))/
               (np.sqrt(pi)*np.cbrt(3*V_d)) + 8*K*C_i*np.sqrt(pi*S_i))/(3*eta)
          
          return P
      

      • testing
        

        x_range = np.logspace(6,9)
        boundary = kinetic_spont_P(V_d=x_range)
        boundary_test(x_param='V_d', y_param='P', x_range=x_range, boundary=boundary,
                      y_offset=-0.5)
        

    • DONE \(\sigma\)
      

      ---------------------------------------------------------------------------
      NameError                                 Traceback (most recent call last)
      <ipython-input-35-594c148ecf9a> in <module>
      ----> 1 def kinetic_spont_sigma(P=P, gamma_m=gamma_m, gamma_d=gamma_d,
            2                   V_d=V_d, S_i=S_i, K=K, C_i=C_i):
            3 
            4     alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
            5     beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
      
      NameError: name 'P' is not defined
      

      • testing
        

        x_range = np.logspace(6,9)
        boundary = kinetic_spont_sigma(V_d=x_range)
        boundary_test(x_param='V_d', y_param='sigma', x_range=x_range, boundary=boundary,
                      y_offset=-0.5)
        

        :END:

    • DONE \(C_{i}\)
      

      from kinetic_spont_ci

      def kinetic_spont_C_i(sigma=sigma, gamma_m=gamma_m, gamma_d=gamma_d,
                        V_d=V_d, S_i=S_i, K=K, P=P):
          
          alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
          eta = (S_i**(3/2))/(6*pi**(1/2))
          
          C_i = (((sigma*alpha + (sigma - gamma_m)*beta +gamma_d*epsilon)*S_i**(3/2))/
               (np.sqrt(pi)*np.cbrt(3*V_d)) + 3*P*eta)/(8*K*np.sqrt(pi*S_i))
          
          return C_i
      

      • testing
        

        x_range = np.logspace(6,9)
        boundary = kinetic_spont_C_i(V_d=x_range)
        boundary_test(x_param='V_d', y_param='C_i', x_range=x_range, boundary=boundary,
                      y_offset=0.5)
        

    • CANCELLED \(\gamma_m\)
      
    • CANCELLED \(\gamma_d\)
      
    • CANCELLED \(S_{i}\)
      
    • DONE \(V_{d}\)
      

      from kinetic_spont_ci

      def kinetic_spont_V_d(sigma=sigma, gamma_m=gamma_m, gamma_d=gamma_d,
                        C_i=C_i, S_i=S_i, K=K, P=P):
          
          alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
          epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
          eta = (S_i**(3/2))/(6*pi**(1/2))
          
          V_d = (1/3)*(((sigma*alpha + (sigma - gamma_m)*beta +gamma_d*epsilon)*S_i**(3/2))/
               ((8*K*C_i*np.sqrt(pi*S_i)-3*P*eta)*np.sqrt(pi)))**3
          
          return V_d
      

      • testing
        

        x_range = np.linspace(0.0002, 0.002)
        boundary = kinetic_spont_V_d(P=x_range)
        boundary_test(x_param='P', y_param='V_d', x_range=x_range, boundary=boundary,
                      y_offset=-5, x_scale='linear')
        

  • F’(0)
    

    def Fprime(h=0, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
               C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i,
               V_d=V_d, P=P):
    
        alpha = np.cbrt((pi*(gamma_m**2-gamma_d**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
        beta = np.cbrt((pi*(gamma_d**2-gamma_m**2)**3)/((2*gamma_d+gamma_m)**2*(gamma_d-gamma_m)**4))
        epsilon = np.cbrt((8*pi)/((2+gamma_m/gamma_d)**2*(1-gamma_m/gamma_d)))
        eta = (S_i**(3/2))/(6*pi**(1/2))
        C = sigma*pi*A_t**2+S_i*(gamma_m-gamma_i)
    
        Fp = ( 3 * P * eta + ( -6 * ( h )**( 2 ) * P * eta + ( -8 * K * ( np.pi )**( 1/2) *
        C_i * ( S_i )**( 1/2 ) + ( ( 1 + -2 * ( h )**( 2 ) ) * (np.pi )**( -1/2 ) * (
        S_i)**( 3/2 ) * ( ( ( 3/2 * h + -1 * ( h )**( 3 ) ) * ( np.pi )**( -1/2 ) * (S_i
        )**( 3/2 ) + 3 * V_d ) )**( -1/3 ) * (sigma*alpha + ( epsilon * gamma_d + beta *
        ( sigma + -1 * gamma_m ) ) ) + 2 * h * (8 * K * np.pi + S_i * (sigma + -1 *
        gamma_m ) ) ) )) )
    
        return Fp
    

  • side-by-side F plots with both phase boundaries
    

    def both_boundary_test(x_param, y_param, x_range, thermo_boundary, kinetic_boundary,
                           params=params, factors=[0.7,1,1.3], y_offset=0.2, x_scale='log', y_scale='linear',
                           variable_df=variable_df, fill=False):
        fparams = params.copy()
        fig, (ax1, ax2)  = plt.subplots(nrows=1, ncols=2, figsize=(12,5))
        
    
        # if x_scale == 'log':
        #     ax1.semilogx(x_range, thermo_boundary, color='black', label='thermodynamic')
        #     ax1.semilogx(x_range, kinetic_boundary, color='blue', label='kinetic')
        # # else:
        #     ax1.plot(x_range, boundary, color='black')
        
        if y_offset > 0:
            fill_range = np.max(kinetic_boundary)
        if y_offset < 0:
            fill_range = np.min(kinetic_boundary)
            
        ax1.loglog(x_range, thermo_boundary, color='orange', label='global stability boundary')
        ax1.loglog(x_range, kinetic_boundary, color='blue', label='local stability boundary')
        if fill == True:
          ax1.fill_between(x_range, y1=thermo_boundary, y2=fill_range*1.3, color='yellow', alpha=0.8, label='globally stable')
          ax1.fill_between(x_range, y1=kinetic_boundary, y2=fill_range*1.3, color='blue', alpha=0.8, label='locally stable')
    
        ax1.legend()
    
        for i, factor in enumerate(factors):
            middle = len(x_range)//2
            x_pick = x_range[middle]*factor
            for j, boundary in enumerate([thermo_boundary, kinetic_boundary]):
              boundary_pick = boundary[middle]+(factor-1)*np.sqrt(thermo_boundary[middle]**2)*y_offset
              fparams[x_param] = x_pick
              fparams[y_param] = boundary_pick
              F = llps_model_recbrt(h=fparams['h_range'], sigma=fparams['sigma'], A_t=fparams['A_t'],
                      S_i=fparams['S_i'], K=fparams['K'], C_i=fparams['C_i'], gamma_d=fparams['gamma_d'],
                      gamma_m=fparams['gamma_m'], gamma_i=fparams['gamma_i'], V_d=fparams['V_d'], P=fparams['P'])
              delta_F = F-F[0]
    
              # ax1.scatter(x_pick, boundary_pick, linewidths=5, color=((j)/2, (j)/2, (j+1)/2, (i+1)/len(factors)))
              # ax2.plot(h_range,delta_F, color=((j)/2, (j)/2, (j+1)/2, (i+1)/len(factors)))
              
              ax1.scatter(x_pick, boundary_pick, linewidths=5)
              ax2.plot(h_range,delta_F)
    
        ax2.hlines(0, -0.05,1.05, color='black')
    
        ax1.set(xlabel=variable_df.loc[x_param]['description']+' ('+variable_df.loc[x_param]['units']+')',
                ylabel=variable_df.loc[y_param]['description']+' ('+variable_df.loc[y_param]['units']+')')
    #            xscale=x_scale)
        ax1.set_xscale(x_scale)
        ax1.set_yscale(y_scale)
        ax2.set(xlabel='h: relative invagination height', ylabel='∆F ($pN \cdot nm$)')
        plt.tight_layout()
    

result plots [8/9]

• plot all outputs vs. invagination height for default parameters
  
  • as individual plots
    

      outputs = llps_model()
      outputs.pop('theta_d')
      # outputs.pop('2÷R_i')
      variables = list(outputs.keys())
    
      for variable in variables:
          plt.figure()
          plot_simulation(output=variable)
          plt.xlabel('h: invagination distance/vesicle diameter')
          plt.savefig(fig_dir+'%s-vs-H_i.png' %(variable))
    

  • as subplots of one figure
    

      outputs = llps_model()
      outputs.pop('theta_d')
      variables = list(outputs.keys())
      nplots = len(variables)
      width = 3
      height = nplots//3+1
      plt.figure(figsize=(width*5,height*4))
    
      for i, variable in enumerate(variables):
          plt.subplot(height, width, i+1)
          plot_simulation(output=variable)
          plt.xlabel('h: invagination distance/vesicle diameter')
    
      plt.tight_layout()
      plt.savefig(fig_dir+'all_outputs.png')
    

• change in all outputs vs. invagination height for default parameters
  
  • as individual plots
    

        outputs = llps_model()
        outputs.pop('theta_d')
        variables = outputs.keys()
    
        for variable in variables:
            plt.figure()
            plot_simulation_delta(output=variable)
            plt.tight_layout()
            plt.savefig(fig_dir+'%s-vs-H_i_delta.png' %(variable))
    

  • as subplots of one figure
    

      outputs = llps_model()
      outputs.pop('theta_d')
      variables = outputs.keys()
      nplots = len(variables)
      width = 3
      height = nplots//3+1
      plt.figure(figsize=(width*5,height*4))
    
      for i, variable in enumerate(variables):
          plt.subplot(height, width, i+1)
          plot_simulation_delta(output=variable)
    
      plt.tight_layout()
      plt.savefig(fig_dir+'all_outputs_delta.png')
    

• normalized geometry terms
  

  plot_params = [
      ['2÷R_i',r'$\frac{1}{R_i}$ $(nm^{-1})$','purple','-'],
      ['V_i','$V_i$ $(nm^3)$','brown','-'],
      ['S_d','$S_d$ $(nm^2)$','green',':'],
      ['S_m','$S_m$ $(nm^2)$','cyan','-'],
      ['S_o','$S_o$ $(nm^2)$','blue','-']
  ]
  
  plt.figure(figsize=(6,6))
  
  for variable, label, color, linestyle in plot_params:
      plot_simulation_normalized(output=variable, H_i=H_i_range,
                      sigma=model_params['sigma'], A_t=model_params['A_t'],
                      S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                      gamma_d=model_params['gamma_d'],
                      gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                      gamma_i=model_params['gamma_1'], V_d=model_params['Vd'],
                      P=model_params['P'], label=label, color=color, linestyle=linestyle,
                      linewidth=4)
  
  plt.legend()
  #plt.ylim(-100, 100)
  plt.ylabel('normalized shape change')
  #plt.title(title)
  plt.tight_layout()
  plt.savefig(fig_dir+'normalized_shapes.png')
  plt.savefig(fig_dir+'normalized_shapes.svg')
  

• change in surface areas
  

  plot_params = [
      ['S_d','$S_d$','green','-'],
      ['S_m','$S_m$','cyan','-']
  ]
  
  plt.figure(figsize=(6,4))
  
  for variable, label, color, linestyle in plot_params:
      plot_simulation_delta(output=variable,
                      sigma=model_params['sigma'], A_t=model_params['A_t'],
                      S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                      gamma_d=model_params['gamma_d'],
                      gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                      gamma_i=model_params['gamma_1'], V_d=model_params['Vd'],
                      P=model_params['P'], label=label, color=color, linestyle=linestyle,
                      linewidth=4)
  plt.ylabel('∆ surface area $(nm^2)$')
  plt.legend()
  plt.tight_layout()
  plt.savefig(fig_dir+'sd_si_vs-H_i_delta.svg')
  plt.savefig(fig_dir+'sd_si_vs-H_i_delta.png')
  
  plt.figure(figsize=(6,4))
  
  for variable, label, color, linestyle in plot_params:
      plot_simulation_delta(x=0.3*H_i_range/vesicle_diameter, H_i=0.3*H_i_range,output=variable,
                      sigma=model_params['sigma'], A_t=model_params['A_t'],
                      S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                      gamma_d=model_params['gamma_d'],
                      gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                      gamma_i=model_params['gamma_1'], V_d=model_params['Vd'],
                      P=model_params['P'], label=label, color=color, linestyle=linestyle,
                      linewidth=4)
  plt.ylabel('∆ surface area $(nm^2)$')
  plt.legend()
  plt.tight_layout()
  plt.savefig(fig_dir+'sd_si_vs-H_i_delta_zoom.svg')
  plt.savefig(fig_dir+'sd_si_vs-H_i_delta_zoom.png')
  
  plot_params = [
      ['surface_tension','$\gamma_d S_d$','green','-'],
      ['wetting_energy','$- \gamma_m S_m$','cyan','-'],
      ['droplet_energy','$F_d$','black','-']
  ]
  
  plt.figure(figsize=(6,4))
  plt.hlines(0,0,1,linestyles=':',colors='gray')
  
  for variable, label, color, linestyle in plot_params:
      plot_simulation_delta(output=variable,
                      sigma=model_params['sigma'], A_t=model_params['A_t'],
                      S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                      gamma_d=model_params['gamma_d'],
                      gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                      gamma_i=model_params['gamma_1'], V_d=model_params['Vd'],
                      P=model_params['P'], label=label, color=color, linestyle=linestyle,
                      linewidth=4)
  plt.ylabel('∆ energy $(pN \cdot nm)$')
  plt.legend()
  plt.tight_layout()
  plt.savefig(fig_dir+'droplet_energies_vs-H_i_delta.svg')
  plt.savefig(fig_dir+'droplet_energies_vs-H_i_delta.png')
  

• each individual energy term output
  
  • set temporary variables
    

      P_temp = 1e-3
      gamma_d_temp = gamma_d*10
      gamma_m_temp = gamma_m*10
      K_temp = 160
    

  • normal droplet
    

    show_geometry(H_i=100)
    plt.savefig(fig_dir+'normal_droplet.png')
    

    • with appropriate signs
      

      droplet only:

      variables = ['surface_tension', 'wetting_energy', 'droplet_energy']
      colors = ['blue', 'green', 'black']
      
      plt.figure(figsize=(6,6))
      
      for variable, color in zip(variables, colors):
          plot_simulation_delta(output=variable, label=variable, P=P_temp, gamma_d=gamma_d_temp,
                                gamma_m=gamma_m_temp, K=K_temp, color=color)
          
      plt.legend()
      #plt.ylim(-100, 100)
      plt.ylabel('∆ energy (pN * nm)')
      #plt.title(title)
      plt.tight_layout()
      plt.savefig(fig_dir+'droplet_energy_contributions.png')
      plt.savefig(fig_dir+'droplet_energy_contributions.svg')
      

      all contributions:

      variables = ['surface_tension', 'wetting_energy', 'droplet_energy',
                   'membrane_tension', 'bending_energy', 'turgor_pressure', 'F']
      colors = ['blue', 'green', 'black',
                'red', 'purple', 'brown', 'orange']
      
      plt.figure(figsize=(6,6))
      
      for variable, color in zip(variables, colors):
          plot_simulation_delta(output=variable, label=variable,
                          sigma=model_params['sigma'], A_t=model_params['A_t'],
                          S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                          gamma_d=model_params['gamma_d'],
                          gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                          gamma_i=model_params['gamma_1'], V_d=model_params['Vd'],
                          P=model_params['P'], color=color)
      
      plt.legend()
      #plt.ylim(-100, 100)
      plt.ylabel('∆ energy (pN * nm)')
      #plt.title(title)
      plt.tight_layout()
      plt.savefig(fig_dir+'all_energy_contributions.png')
      plt.savefig(fig_dir+'all_energy_contributions.svg')
      

      variables = ['membrane_tension', 'bending_energy',
                   'surface_tension', 'wetting_energy',
                   'turgor_pressure']
      labels = ['membrane tension', 'bending energy',
                   'droplet surface tension', 'wetting energy',
                   'turgor pressure']
      colors = ['red', 'purple', 'blue', 'green', 'brown']
      
      plt.figure(figsize=(6,6))
      
      for variable, label, color in zip(variables, labels, colors):
          plot_simulation_delta(output=variable, label=label, P=1.1e-3, gamma_d=0.06,
                                gamma_m=0.02, K=50, color=color)
      
      plt.legend()
      #plt.ylim(-100, 100)
      plt.ylabel('∆ free energy') # (pN * nm)
      #plt.title(title)
      plt.tight_layout()
      plt.savefig(fig_dir+'individual_energy_contributions.png')
      plt.savefig(fig_dir+'individual_energy_contributions.svg')
      

      effect of droplet:

      variables = ['non_droplet_energy', 'droplet_energy', 'F']
      labels = ['non-condensate energy contributions',
                'condensate energy contributions',
                'total energy']
      colors = ['blue', 'green', 'black']
      
      plt.figure(figsize=(6,5))
      
      for variable, label, color in zip(variables, labels, colors):
          plot_simulation_delta(output=variable, label=label, P=1.1e-3, gamma_d=0.06,
                                gamma_m=0.02, K=50, color=color)
      
      plt.legend()
      #plt.ylim(-100, 100)
      plt.ylabel('∆ free energy ($pN \cdot nm$)') # (pN * nm)
      #plt.title(title)
      plt.tight_layout()
      plt.savefig(fig_dir+'condensate_energy_contribution.png')
      plt.savefig(fig_dir+'condensate_energy_contribution.svg')
      

    • flipped sign of negative energy contributions
      

      variables = ['membrane_tension', 'bending_energy_flipped','surface_tension',
                    'wetting_energy_flipped', 'turgor_pressure']
      
      plt.figure(figsize=(10,10))
      
      for variable in variables:
          plot_simulation_delta(output=variable, label=variable, P=P_temp, gamma_d=gamma_d_temp,
                                gamma_m=gamma_m_temp, K=K_temp)
      plt.legend()
      #plt.ylim(-100, 100)
      plt.ylabel('∆ energy (pN * nm)')
      #plt.title(title)
      plt.tight_layout()
      plt.savefig(fig_dir+'energy_contributions_flipped.png')
      

  • huge droplet
    

    show_geometry(H_i=100, V_d=V_d*10**4,
                  left=-3600, right=3600, bottom=-2300)
    plt.savefig(fig_dir+'bigdrop.png')
    

    • with appropriate signs
      

      variables = ['surface_tension', 'wetting_energy', 'droplet_energy']
      colors = ['blue', 'green', 'black']
      
      plt.figure(figsize=(6,6))
      
      for variable, color in zip(variables, colors):
          plot_simulation_delta(output=variable, label=variable, P=P_temp, gamma_d=gamma_d_temp,
                                gamma_m=gamma_m_temp, K=K_temp, V_d=V_d*10**4, color=color)
      plt.legend()
      #plt.ylim(-100, 100)
      plt.ylabel('∆ energy (pN * nm)')
      #plt.title(title)
      plt.tight_layout()
      plt.savefig(fig_dir+'droplet_energy_contributions_bigdrop.png')
      plt.savefig(fig_dir+'droplet_energy_contributions_bigdrop.svg')
      

      variables = ['non_droplet_energy', 'droplet_energy', 'F']
      labels = ['non-condensate energy contributions',
                'condensate energy contributions',
                'total energy']
      colors = ['blue', 'green', 'black']
      
      plt.figure(figsize=(6,5))
      
      for variable, label, color in zip(variables, labels, colors):
          plot_simulation_delta(output=variable, label=label, P=1.1e-3, gamma_d=0.06,
                                gamma_m=0.02, K=50, V_d=V_d*10, color=color)
      
      plt.legend()
      #plt.ylim(-100, 100)
      plt.ylabel('∆ free energy ($pN \cdot nm$)') #pN*nm
      #plt.title(title)
      plt.tight_layout()
      plt.savefig(fig_dir+'condensate_energy_contribution_bigv.png')
      plt.savefig(fig_dir+'condensate_energy_contribution_bigv.svg')
      

    • flipped sign of negative energy contributions
      

      variables = ['membrane_tension', 'bending_energy_flipped','surface_tension',
                    'wetting_energy_flipped', 'turgor_pressure']
      
      plt.figure(figsize=(10,10))
      
      for variable in variables:
          plot_simulation_delta(output=variable, label=variable, P=P_temp, gamma_d=gamma_d_temp,
                                gamma_m=gamma_m_temp, K=K_temp, V_d=V_d*10**4)
      plt.legend()
      #plt.ylim(-100, 100)
      plt.ylabel('∆ energy (pN * nm)')
      #plt.title(title)
      plt.tight_layout()
      plt.savefig(fig_dir+'energy_contributions_bigdrop_flipped.png')
      

• Plot F for a series of values of gamma_m
  

  What are the parameter values? I assume that gamma_m is negative. Plot F for a series of values of gamma_m. Based on the plots for S_m, S_d, and S_0, F should end up with a minimum at H/R=1 if gamma_m is large enough. It might also happen if gamma_d is small enough. It’s hard to say because I can’t easily compare the magnitude of changes in S_d, S_i, and S_0.

  • absolute
    

         gamma_m_range = [0,50,100,150,199]
    
         for i in gamma_m_range:
             plot_simulation(gamma_m=i, label=i/gamma_d)
    
         plt.legend(title='gamma_m/gamma_d')
         plt.xlabel('invagination distance/vesicle diameter')
         plt.tight_layout()
         plt.savefig(fig_dir+'f_scan_gamma_m.png')
    

  • normalized
    

         gamma_m_range = [0,50,100,150,199]
    
         for i in gamma_m_range:
             plot_simulation_normalized(gamma_m=i, label=i/gamma_d)
    
         plt.legend(title='gamma_m/gamma_d')
         plt.xlabel('invagination distance/vesicle diameter')
         plt.tight_layout()
         plt.savefig(fig_dir+'f_scan_gamma_m.png')
    

    :end:

  • delta
    

       gamma_m_range = np.array( [0.01,0.25,0.5,0.75] )*gamma_d
    
       for i in gamma_m_range:
           plot_simulation_delta(gamma_m=i, label=i/gamma_d)
    
       plt.legend(title='$\gamma_{m}/\gamma_{d}$')
       plt.tight_layout()
       plt.savefig(fig_dir+'f_scan_gamma_m_delta.png')
    

  • delta and geometry side by side
    

         gamma_m_range = [0,50,100,150,180]
    
         for i in gamma_m_range:
             plot_simulation_geometry(gamma_m=i, label=i/gamma_d)
             plt.legend(title='$\gamma_{m}/\gamma_{d}$')
             plt.tight_layout()
             plt.savefig(fig_dir+'f_scan_gamma_m_%s.png' %(i))
    

• F over series of turgor pressure
  

       pressure_range = np.linspace(10**-2,1,5)
  
       for i in pressure_range:
           plot_simulation_delta(P=i, label=i)
  
       plt.legend(title='turgor pressure ($N/nm^{2}$)')
       plt.tight_layout()
       plt.savefig(fig_dir+'f_scan_turgor_delta.png')
  

• effect of droplet volume on energy barrier
  

  #volume_range = np.linspace(2/3*pi*200**3, 2/3*pi*2000**3,1000)
  volume_range = np.logspace(7, 11, 100)
  fc_prime_0 = Fprime(h=0, sigma=0, A_t=A_t, S_i=S_i, K=0,
             C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i,
             V_d=volume_range, P=0)
  
  plt.figure(figsize=(7,6))
  plt.plot(volume_range,fc_prime_0)
  plt.xscale('log')
  #plt.yscale('log')
  plt.xlabel('condensate volume ($nm^3$)')
  plt.ylabel(r'$\frac{dF_c}{dh} (h=0)$'+'\n (slope of condensate free energy at h=0)')
  plt.tight_layout()
  plt.savefig(fig_dir+'dfc_dh_0_vdscan.svg')
  plt.savefig(fig_dir+'dfc_dh_0_vdscan.png')
  

  #volume_range = np.linspace(2/3*pi*200**3, 2/3*pi*2000**3,1000)
  volume_range = np.logspace(7, 11, 100)
  full_fc_prime_0 = Fprime(h=0, sigma=0, A_t=A_t, S_i=S_i, K=0,
             C_i=C_i, gamma_d=gamma_d, gamma_m=gamma_m, gamma_i=gamma_i,
             V_d=volume_range, P=0)
  tension_fc_prime_0 = Fprime(h=0, sigma=0, A_t=A_t, S_i=S_i, K=0,
             C_i=C_i, gamma_d=gamma_d, gamma_m=10e-20, gamma_i=10e-20,
             V_d=volume_range, P=0)
  wetting_fc_prime_0 = Fprime(h=0, sigma=0, A_t=A_t, S_i=S_i, K=0,
             C_i=C_i, gamma_d=10e-20, gamma_m=gamma_m, gamma_i=10e-20,
             V_d=volume_range, P=0)
  
  plt.figure(figsize=(7,6))
  plt.plot(volume_range,tension_fc_prime_0,label='surface tension')
  plt.plot(volume_range,wetting_fc_prime_0,label='wetting energy')
  plt.plot(volume_range,full_fc_prime_0,label='sum')
  plt.xscale('log')
  #plt.yscale('log')
  plt.legend()
  plt.xlabel('condensate volume ($nm^3$)')
  plt.ylabel(r'$\frac{dF_c}{dh} (h=0)$'+'\n (slope of condensate free energy at h=0)')
  plt.tight_layout()
  plt.savefig(fig_dir+'dfcs_dh_0_vdscan.svg')
  plt.savefig(fig_dir+'dfcs_dh_0_vdscan.png')
  

  volume_range = np.logspace(7, 12, 100)
  dfmaxlist = []
  
  for v in volume_range:
      f = llps_model_recbrt(h=h_range, sigma=sigma, A_t=A_t, S_i=S_i, K=K,
                  C_i=C_i, gamma_d=0.5, gamma_m=0.2, gamma_i=0.02,
                  V_d=v, P=0.005)
      dfmax = f.max()-f[0]
      dfmaxlist.append(dfmax)
  
  plt.figure(figsize=(7,6))
  plt.plot(volume_range,dfmaxlist)
  plt.xscale('log')
  #plt.yscale('log')
  plt.xlabel('condensate volume ($nm^3$)')
  plt.ylabel(r'$\Delta F_{max}$')
  plt.tight_layout()
  plt.savefig(fig_dir+'dfmax_vdscan.svg')
  plt.savefig(fig_dir+'dfmax_vdscan.png')
  

  • plot energy derivatives vs. volume
    

    droplet only

    V_d_range = np.linspace(5E6,5E7)
    
    plt.plot(V_d_range,Fprime(h=0,
                        sigma=0, A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=0, C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'],
                        gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=V_d_range,
                              P=0), label=r'$\frac{d F_{d}}{d h}$', color='black',# linestyle=linestyle,
                        linewidth=4)
    
    # plt.plot(V_d_range,Fprime(h=0,
    #                     sigma=model_params['sigma'], A_t=model_params['A_t'],
    #                     S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
    #                     gamma_d=model_params['gamma_d'],
    #                     gamma_m=model_params['gamma_3']-model_params['gamma_2'],
    #                     gamma_i=model_params['gamma_1'], V_d=V_d_range,
    #                     P=model_params['P']), label='F all', color='orange', linestyle=linestyle,
    #                     linewidth=4)
    
    plt.plot(V_d_range,-Fprime(h=0,
                        sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'],
                        gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=V_d_range,
                        P=model_params['P']) +
             Fprime(h=0,
                        sigma=0, A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=0, C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'],
                        gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=V_d_range,
                              P=0), label=r'$-\frac{d F_{e}}{d h}$', color='orange',# linestyle=linestyle,
                        linewidth=4)
    
    plt.legend()
    plt.xlabel('$V_d$: droplet volume')
    plt.ylabel('slope of free energy at $h=0$')
    plt.tight_layout()
    plt.savefig(fig_dir+'dfdh_vdscan.svg')
    plt.savefig(fig_dir+'dfdh_vdscan.png')
    

  • animation
    

    animate_simulation_geometry(
                        sigma=model_params['sigma'], A_t=model_params['A_t']*1000,
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'],
                        gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d_range=V_d_range,
                        P=model_params['P'],
                    geometry_top=20, geometry_bottom=-270, geometry_left=-400, geometry_right=400,
                    plot_top=2300, plot_bottom=-1900, plot_left=-0.1, plot_right=0.1
                        )
    

• spontaneity analysis
  

  there’s no solution for kinetic \(\sigma\)

  • DONE \(V_{d}\) x-axis
    

    boundaries = [(thermo_spont_K,kinetic_spont_K, 'K', 2),
                  (thermo_spont_P,kinetic_spont_P, 'P', -2),
                  #(thermo_spont_sigma,kinetic_spont_sigma, 'sigma'),
                  (thermo_spont_C_i,kinetic_spont_C_i, 'C_i', 2)]
    x_range = np.logspace(7,9,base=10)
    for thermo_boundaryy, kinetic_boundaryy, y_param, y_offset in boundaries:
        thermo_boundary = thermo_boundaryy(V_d=x_range)
        kinetic_boundary = kinetic_boundaryy(V_d=x_range)
        both_boundary_test(x_param='V_d', y_param=y_param, x_range=x_range,
                        thermo_boundary=thermo_boundary, kinetic_boundary=kinetic_boundary,
                        y_offset=y_offset/2, factors=[0.9,1,1.1])
        plt.savefig(fig_dir+y_param+'V_d_boundary_check.png')
    

    show_geometry(H_i=100, V_d=10**(7))
    plt.savefig(fig_dir+'10e7_geometry.png')
    

    show_geometry(H_i=100, V_d=10**(9), left=-1200, right=1200, bottom=-600)
    plt.savefig(fig_dir+'10e9_geometry.png')
    

  • DONE \(\gamma_{d}\) x-axis
    

    boundaries = [(thermo_spont_K,kinetic_spont_K, 'K', -4),
                  (thermo_spont_P,kinetic_spont_P, 'P', 2),
                  #(thermo_spont_sigma,kinetic_spont_sigma, 'sigma'),
                  (thermo_spont_C_i,kinetic_spont_C_i, 'C_i', 2)]
    x_range = np.logspace(-1,0,base=10)
    for thermo_boundaryy, kinetic_boundaryy, y_param, y_offset in boundaries:
        thermo_boundary = thermo_boundaryy(gamma_d=x_range)
        kinetic_boundary = kinetic_boundaryy(gamma_d=x_range)
        both_boundary_test(x_param='gamma_d', y_param=y_param, x_range=x_range,
                        thermo_boundary=thermo_boundary, kinetic_boundary=kinetic_boundary,
                        y_offset=y_offset/2, factors=[0.95,1,1.05])
        plt.savefig(fig_dir+y_param+'gamma_d_boundary_check.png')
        plt.savefig(fig_dir+y_param+'gamma_d_boundary_check.svg')
    

  • DONE \(\gamma_{m}\) x-axis
    

    boundaries = [(thermo_spont_K,kinetic_spont_K, 'K', 2),
                  (thermo_spont_P,kinetic_spont_P, 'P', -2),
                  #(thermo_spont_sigma,kinetic_spont_sigma, 'sigma'),
                  (thermo_spont_C_i,kinetic_spont_C_i, 'C_i', 2)]
    x_range = np.logspace(-3,0,base=10)
    for thermo_boundaryy, kinetic_boundaryy, y_param, y_offset in boundaries:
        thermo_boundary = thermo_boundaryy(gamma_m=x_range)
        kinetic_boundary = kinetic_boundaryy(gamma_m=x_range)
        both_boundary_test(x_param='gamma_m', y_param=y_param, x_range=x_range,
                        thermo_boundary=thermo_boundary, kinetic_boundary=kinetic_boundary,
                        y_offset=y_offset*2, factors=[0.95,1,1.05])
        plt.savefig(fig_dir+y_param+'gamma_m_boundary_check.png')
    

  • TODO \(\gamma_{m}\) /\(\gamma_{d}\) x-axis
    

    boundaries = [(thermo_spont_K,kinetic_spont_K, 'K', 2),
                  (thermo_spont_P,kinetic_spont_P, 'P', -2),
                  #(thermo_spont_sigma,kinetic_spont_sigma, 'sigma'),
                  (thermo_spont_C_i,kinetic_spont_C_i, 'C_i', 2)]
    x_range = np.logspace(-2,0.0001,base=10)*gamma_d
    for thermo_boundaryy, kinetic_boundaryy, y_param, y_offset in boundaries:
        thermo_boundary = thermo_boundaryy(gamma_m=x_range)
        kinetic_boundary = kinetic_boundaryy(gamma_m=x_range)
        both_boundary_test(x_param='gamma_m', y_param=y_param, x_range=x_range,
                        thermo_boundary=thermo_boundary, kinetic_boundary=kinetic_boundary,
                        y_offset=y_offset/2, factors=[0.8,1,1.2])
        plt.savefig(fig_dir+y_param+'gamma_m_relative_boundary_check.png')
    

  • DONE \(\sigma\) x-axis
    

    boundaries = [(thermo_spont_K,kinetic_spont_K, 'K', -2),
                  (thermo_spont_P,kinetic_spont_P, 'P', 2),
                  (thermo_spont_C_i,kinetic_spont_C_i, 'C_i', -2)]
    x_range = np.logspace(-5,-1,base=10)
    for thermo_boundaryy, kinetic_boundaryy, y_param, y_offset in boundaries:
        thermo_boundary = thermo_boundaryy(sigma=x_range)
        kinetic_boundary = kinetic_boundaryy(sigma=x_range)
        both_boundary_test(x_param='sigma', y_param=y_param, x_range=x_range,
                        thermo_boundary=thermo_boundary, kinetic_boundary=kinetic_boundary,
                        y_offset=y_offset/5, factors=[0.8,1,1.2])
        plt.savefig(fig_dir+y_param+'sigma_boundary_check.png')
    

  • DONE \(S_{i}\) x-axis
    

    boundaries = [(thermo_spont_K,kinetic_spont_K, 'K', 2),
                  (thermo_spont_P,kinetic_spont_P, 'P', -2),
                  (thermo_spont_C_i,kinetic_spont_C_i, 'C_i', 2)]
    x_range = np.logspace(3.5,4.5,base=10)
    for thermo_boundaryy, kinetic_boundaryy, y_param, y_offset in boundaries:
        thermo_boundary = thermo_boundaryy(S_i=x_range)
        kinetic_boundary = kinetic_boundaryy(S_i=x_range)
        both_boundary_test(x_param='S_i', y_param=y_param, x_range=x_range,
                        thermo_boundary=thermo_boundary, kinetic_boundary=kinetic_boundary,
                        y_offset=y_offset*50, factors=[0.999,1,1.001])
        plt.savefig(fig_dir+y_param+'S_i_boundary_check.png')
    

  • DONE \(P\) x-axis
    

    boundaries = [(thermo_spont_K,kinetic_spont_K, 'K', 2),
                  #(thermo_spont_P,kinetic_spont_P, 'P', 2),
                  (thermo_spont_C_i,kinetic_spont_C_i, 'C_i', 2)]
    x_range = np.logspace(-3,0,base=10)
    for thermo_boundaryy, kinetic_boundaryy, y_param, y_offset in boundaries:
        thermo_boundary = thermo_boundaryy(P=x_range)
        kinetic_boundary = kinetic_boundaryy(P=x_range)
        both_boundary_test(x_param='P', y_param=y_param, x_range=x_range,
                        thermo_boundary=thermo_boundary, kinetic_boundary=kinetic_boundary,
                        y_offset=y_offset*150, factors=[0.999,1,1.001])
        plt.savefig(fig_dir+y_param+'P_boundary_check.png')
    

  • DONE \(C_{i}\) x-axis
    

    boundaries = [(thermo_spont_K,kinetic_spont_K, 'K', 2),
                  (thermo_spont_P,kinetic_spont_P, 'P', -2)]
    x_range = np.logspace(-4,-1,base=10)
    for thermo_boundaryy, kinetic_boundaryy, y_param, y_offset in boundaries:
        thermo_boundary = thermo_boundaryy(C_i=x_range)
        kinetic_boundary = kinetic_boundaryy(C_i=x_range)
        both_boundary_test(x_param='C_i', y_param=y_param, x_range=x_range,
                        thermo_boundary=thermo_boundary, kinetic_boundary=kinetic_boundary,
                        y_offset=y_offset/5, factors=[0.8,1,1.2])
        plt.savefig(fig_dir+y_param+'C_i_boundary_check.png')
    

  • DONE \(K\) x-axis
    

    boundaries = [(thermo_spont_C_i,kinetic_spont_C_i, 'C_i', 2),
                  (thermo_spont_P,kinetic_spont_P, 'P', -2)]
    x_range = np.logspace(1,4,base=10)
    for thermo_boundaryy, kinetic_boundaryy, y_param, y_offset in boundaries:
        thermo_boundary = thermo_boundaryy(K=x_range)
        kinetic_boundary = kinetic_boundaryy(K=x_range)
        both_boundary_test(x_param='K', y_param=y_param, x_range=x_range,
                        thermo_boundary=thermo_boundary, kinetic_boundary=kinetic_boundary,
                        y_offset=y_offset/1, factors=[0.9,1,1.1])
        plt.savefig(fig_dir+y_param+'K_boundary_check.png')
    

• plots for 20200915 lab meeting drafting
  
  • TODO each individual energy term output
    

        variables = ['membrane_tension',
                      'bending_energy','surface_tension',
                      'wetting_energy','turgor_pressure']
    
        for variable in variables:
            plt.figure()
            plot_simulation_delta(output=variable)
            plt.tight_layout()
            plt.savefig(fig_dir+'%s-vs-H_i_delta.png' %(variable))
    

  • membrane components without droplet
    
    • effect of membrane tension alone
      

       plot_simulation_geometry(V_d=0, K=0, gamma_m=0, gamma_i=0, P=0)
      

    • effect of bending energy alone
      

       plot_simulation_geometry(K=K, V_d=0, sigma=0, gamma_m=0, gamma_i=0, P=0)
      

      • TODO fix this
        

        why is this different from plotting membrane bending energy term alone?

         R_i_range = np.linspace(0.1,np.pi,100)
         bending_energy = (K/2)*S_i*((2/R_i_range)-C_i)**2
         plt.plot(R_i_range, bending_energy)
        

  • droplet effects
    
    • surface tension
      

       plot_simulation_geometry(sigma=0, K=0, gamma_m=0, gamma_i=0, P=0)
      

    • surface tension + wetting energy
      

      can’t have wetting energy alone

       plot_simulation_geometry(sigma=0, K=0, P=0)
      

  • TODO droplet+membrane effects
    
  • TODO turgor pressure effects
    
• complementing membrane mechanics model
  
  • kinetically spontaneous
    

    N=10
    #coat_area = np.pi*100**2
    #coat_radius = np.sqrt(coat_area/(4*np.pi))
    #total_radius = np.sqrt(coat_area/np.pi)*3
    #Vd = 1/3*np.pi*200**3
    
    coat_radius = 50
    #coat_area = ((3*N)/(3*N+2))*4*np.pi*coat_radius**2
    coat_area = 4*np.pi*coat_radius**2
    total_radius = np.sqrt(coat_area/np.pi)*3
    Vd = 1/2*np.pi*200**3
    
    model_params = {
        'time':1000000,
        'dt':1,
        'print_interval':10000,
        'write_interval':1000,
        'shape':'initiation',
        'P':0.0001,
        'sigma':0.01,
        'gamma_1':0,
        'gamma_2':0.002,
        'gamma_3':0.02,
        'gamma_d':0.04,
        'h1':1E-2,
        'h2':1E-2,
        'h3':1E-2,
        'H1':1/coat_radius,
        'K1':50,
        'H3':0,
        'K3':1,
        'polar_type':'integral',
        'G':1,
        'cell_wall':True,
        'N':N,
        'A_t':total_radius,
        'A_c':coat_area,
        'Vd':Vd,
        'emap':False,
        'tip_pull':0
    }
    
    

    plot_simulation_geometry(x=H_i_range/vesicle_diameter, output='F',
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         fig_width=15, fig_height=5, label=None, variable_df=variable_df,
                        num_plots=2, top=20, bottom=-250, left=-400, right=400)
    

    variables = ['surface_tension', 'wetting_energy', 'droplet_energy',
                 'membrane_tension', 'bending_energy', 'turgor_pressure', 'F']
    colors = ['blue', 'green', 'black',
              'red', 'purple', 'brown', 'orange']
    
    plt.figure(figsize=(6,6))
    
    for variable, color in zip(variables, colors):
        plot_simulation_delta(output=variable, label=variable, color=color,
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         variable_df=variable_df)
    
    plt.legend()
    plt.hlines(0,-0.1,1.1,linestyles='dotted',colors='black')
    #plt.ylim(-100, 100)
    plt.ylabel('∆ energy (pN * nm)')
    #plt.title(title)
    plt.tight_layout()
    plt.savefig(fig_dir+'energy_contributions_complementary.png')
    plt.savefig(fig_dir+'energy_contributions_complementary.svg')
    

    plot_simulation_delta(output='F', label='F', color='orange',
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         variable_df=variable_df)
    plt.legend()
    #plt.hlines(0,-0.1,1.1,linestyles='dotted',colors='black')
    #plt.ylim(-100, 100)
    plt.ylabel('∆ energy (pN * nm)')
    #plt.title(title)
    plt.tight_layout()
    plt.savefig(fig_dir+'f_complementary.png')
    plt.savefig(fig_dir+'f_complementary.svg')
    

  • only thermodynamically spontaneous
    

    N=10
    #coat_area = np.pi*100**2
    #coat_radius = np.sqrt(coat_area/(4*np.pi))
    #total_radius = np.sqrt(coat_area/np.pi)*3
    #Vd = 1/3*np.pi*200**3
    
    coat_radius = 50
    #coat_area = ((3*N)/(3*N+2))*4*np.pi*coat_radius**2
    coat_area = 4*np.pi*coat_radius**2
    total_radius = np.sqrt(coat_area/np.pi)*3
    Vd = 1/2*np.pi*200**3
    
    model_params = {
        'time':1000000,
        'dt':1,
        'print_interval':10000,
        'write_interval':1000,
        'shape':'initiation',
        'P':0.0001,
        'sigma':0.01,
        'gamma_1':0,
        'gamma_2':0.002,
        'gamma_3':0.02,
        'gamma_d':0.06,
        'h1':1E-2,
        'h2':1E-2,
        'h3':1E-2,
        'H1':1/coat_radius,
        'K1':10,
        'H3':0,
        'K3':1,
        'polar_type':'integral',
        'G':1,
        'cell_wall':True,
        'N':N,
        'A_t':total_radius,
        'A_c':coat_area,
        'Vd':Vd,
        'emap':False,
        'tip_pull':0
    }
    
    

    plot_simulation_geometry(x=H_i_range/vesicle_diameter, output='F',
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         fig_width=15, fig_height=5, label=None, variable_df=variable_df,
                        num_plots=2, top=20, bottom=-250, left=-400, right=400)
    

    variables = ['surface_tension', 'wetting_energy', 'droplet_energy',
                 'membrane_tension', 'bending_energy', 'turgor_pressure', 'F']
    colors = ['blue', 'green', 'black',
              'red', 'purple', 'brown', 'orange']
    
    plt.figure(figsize=(6,6))
    
    for variable, color in zip(variables, colors):
        plot_simulation_delta(output=variable, label=variable, color=color,
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         variable_df=variable_df)
    
    plt.legend()
    plt.hlines(0,-0.1,1.1,linestyles='dotted',colors='black')
    #plt.ylim(-100, 100)
    plt.ylabel('∆ energy (pN * nm)')
    #plt.title(title)
    plt.tight_layout()
    plt.savefig(fig_dir+'energy_contributions_complementary_barrier.png')
    plt.savefig(fig_dir+'energy_contributions_complementary_barrier.svg')
    

    plot_simulation_delta(output='F', label='F', color='orange',
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         variable_df=variable_df)
    plt.legend()
    #plt.hlines(0,-0.1,1.1,linestyles='dotted',colors='black')
    #plt.ylim(-100, 100)
    plt.ylabel('∆ energy (pN * nm)')
    #plt.title(title)
    plt.tight_layout()
    plt.savefig(fig_dir+'f_complementary_barrier.png')
    plt.savefig(fig_dir+'f_complementary_barrier.svg')
    

  • not spontaneous
    

    N=10
    #coat_area = np.pi*100**2
    #coat_radius = np.sqrt(coat_area/(4*np.pi))
    #total_radius = np.sqrt(coat_area/np.pi)*3
    #Vd = 1/3*np.pi*200**3
    
    coat_radius = 50
    #coat_area = ((3*N)/(3*N+2))*4*np.pi*coat_radius**2
    coat_area = 4*np.pi*coat_radius**2
    total_radius = np.sqrt(coat_area/np.pi)*3
    Vd = 1/2*np.pi*200**3
    
    model_params = {
        'time':1000000,
        'dt':1,
        'print_interval':10000,
        'write_interval':1000,
        'shape':'initiation',
        'P':0.0001,
        'sigma':0.01,
        'gamma_1':0,
        'gamma_2':0.002,
        'gamma_3':0.02,
        'gamma_d':0.06,
        'h1':1E-2,
        'h2':1E-2,
        'h3':1E-2,
        'H1':1/coat_radius,
        'K1':1,
        'H3':0,
        'K3':1,
        'polar_type':'integral',
        'G':1,
        'cell_wall':True,
        'N':N,
        'A_t':total_radius,
        'A_c':coat_area,
        'Vd':Vd,
        'emap':False,
        'tip_pull':0
    }
    
    

    plot_simulation_geometry(x=H_i_range/vesicle_diameter, output='F',
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         fig_width=15, fig_height=5, label=None, variable_df=variable_df,
                        num_plots=2, top=20, bottom=-250, left=-400, right=400)
    

    variables = ['surface_tension', 'wetting_energy', 'droplet_energy',
                 'membrane_tension', 'bending_energy', 'turgor_pressure', 'F']
    colors = ['blue', 'green', 'black',
              'red', 'purple', 'brown', 'orange']
    
    plt.figure(figsize=(6,6))
    
    for variable, color in zip(variables, colors):
        plot_simulation_delta(output=variable, label=variable, color=color,
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         variable_df=variable_df)
    
    plt.legend()
    plt.hlines(0,-0.1,1.1,linestyles='dotted',colors='black')
    #plt.ylim(-100, 100)
    plt.ylabel('∆ energy (pN * nm)')
    #plt.title(title)
    plt.tight_layout()
    plt.savefig(fig_dir+'energy_contributions_complementary_notspon.png')
    plt.savefig(fig_dir+'energy_contributions_complementary_notspon.svg')
    

    plot_simulation_delta(output='F', label='F', color='orange',
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         variable_df=variable_df)
    plt.legend()
    #plt.hlines(0,-0.1,1.1,linestyles='dotted',colors='black')
    #plt.ylim(-100, 100)
    plt.ylabel('∆ energy (pN * nm)')
    #plt.title(title)
    plt.tight_layout()
    plt.savefig(fig_dir+'f_complementary_notspon.png')
    plt.savefig(fig_dir+'f_complementary_notspon.svg')
    

  • kinetically spontaneous flat
    

    N=10
    #coat_area = np.pi*100**2
    #coat_radius = np.sqrt(coat_area/(4*np.pi))
    #total_radius = np.sqrt(coat_area/np.pi)*3
    #Vd = 1/3*np.pi*200**3
    
    coat_radius = 50
    #coat_area = ((3*N)/(3*N+2))*4*np.pi*coat_radius**2
    coat_area = 4*np.pi*coat_radius**2
    total_radius = np.sqrt(coat_area/np.pi)*3
    Vd = 1/2*np.pi*200**3
    
    model_params = {
        'time':1000000,
        'dt':1,
        'print_interval':10000,
        'write_interval':1000,
        'shape':'initiation',
        'P':0.0001,
        'sigma':0.1,
        'gamma_1':0,
        'gamma_2':0.002,
        'gamma_3':0.02,
        'gamma_d':0.06,
        'h1':1E-2,
        'h2':1E-2,
        'h3':1E-2,
        'H1':1/coat_radius,
        'K1':1,
        'H3':0,
        'K3':1,
        'polar_type':'integral',
        'G':1,
        'cell_wall':True,
        'N':N,
        'A_t':total_radius,
        'A_c':coat_area,
        'Vd':Vd,
        'emap':False,
        'tip_pull':0
    }
    
    

    plot_simulation_geometry(x=H_i_range/vesicle_diameter, output='F',
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         fig_width=15, fig_height=5, label=None, variable_df=variable_df,
                        num_plots=2, top=20, bottom=-250, left=-400, right=400)
    

    variables = ['surface_tension', 'wetting_energy', 'droplet_energy',
                 'membrane_tension', 'bending_energy', 'turgor_pressure', 'F']
    colors = ['blue', 'green', 'black',
              'red', 'purple', 'brown', 'orange']
    
    plt.figure(figsize=(6,6))
    
    for variable, color in zip(variables, colors):
        plot_simulation_delta(output=variable, label=variable, color=color,
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         variable_df=variable_df)
    
    plt.legend()
    plt.hlines(0,-0.1,1.1,linestyles='dotted',colors='black')
    #plt.ylim(-100, 100)
    plt.ylabel('∆ energy (pN * nm)')
    #plt.title(title)
    plt.tight_layout()
    plt.savefig(fig_dir+'energy_contributions_complementary_sponflat.png')
    plt.savefig(fig_dir+'energy_contributions_complementary_sponflat.svg')
    

    plot_simulation_delta(output='F', label='F', color='orange',
                         H_i=H_i_range, sigma=model_params['sigma'], A_t=model_params['A_t'],
                        S_i=model_params['A_c'], K=model_params['K1'], C_i=model_params['H1'],
                        gamma_d=model_params['gamma_d'], gamma_m=model_params['gamma_3']-model_params['gamma_2'],
                        gamma_i=model_params['gamma_1'], V_d=model_params['Vd'], P=model_params['P'],
                         variable_df=variable_df)
    plt.legend()
    #plt.hlines(0,-0.1,1.1,linestyles='dotted',colors='black')
    #plt.ylim(-100, 100)
    plt.ylabel('∆ energy (pN * nm)')
    #plt.title(title)
    plt.tight_layout()
    plt.savefig(fig_dir+'f_complementary_sponflat.png')
    plt.savefig(fig_dir+'f_complementary_sponflat.svg')